Section 4 Scenario Analysis: WCS Priority Landscapes

# turn off the s2 processing 
## https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data
sf::sf_use_s2(FALSE)

See the USFS Wildfire Crisis Strategy (WCS) document for background information.

4.1 Scenarios Considered

There are three scenarios considered in this analysis of the constraints to mechanical forest health and fuel reduction treatments:

# mechanical mgmt allowed if:
n_sc_temp <- 3
data.frame(
  scenario = 1:3
  , cover = rep("Forest or Shrubland",n_sc_temp)
  , protected = rep("Not Protected",n_sc_temp)
  , slope = c("<40%","<60%","<60%")
  , road = c("<1,000ft","<2,000ft","<2,000ft")
  , riparian = c(">100ft",">100ft",">50ft")
  , admin = c("No Designation","No Designation","Any Designation")
) %>% 
  tidyr::pivot_longer(
    cols = -c(scenario)
  ) %>% 
  tidyr::pivot_wider(
    names_from = scenario
    , values_from = value
    , names_prefix = "scenario_"
  ) %>% 
  dplyr::mutate(
    name = factor(
      name
      , levels = c(
        "cover"
        , "protected"
        , "slope"
        , "road"
        , "riparian"
        , "admin"
      )
      , labels = c(
        "NLCD Cover Type"
        , "Protected or<br>IRA Status"
        , "Slope"
        , "Distance to<br>Nearest Road"
        , "Riparian<br>Buffer"
        , "Administrative<br>Designation"
      )
      , ordered = T
    )
  ) %>% 
  dplyr::arrange(name) %>% 
  dplyr::mutate(
    lvl = paste0("L",dplyr::row_number()-1,":")
  ) %>% 
  dplyr::relocate(lvl) %>% 
# make table
kableExtra::kable(
    caption = "Scenarios Considered for Determining Mechanical Treatment Operability<br>Mechanical treatment allowed if:"
    , col.names = c(
      ""
      , ""
      , "Scenario 1<br>Most Constrained"
      , "Scenario 2<br>Moderate"
      , "Scenario 3<br>Least Constrained"
    )
    , escape = F
  ) %>%  
  kable_classic(full_width=T) %>% 
  kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE)
Table 4.1: Scenarios Considered for Determining Mechanical Treatment Operability
Mechanical treatment allowed if:
Scenario 1
Most Constrained
Scenario 2
Moderate
Scenario 3
Least Constrained
L0: NLCD Cover Type Forest or Shrubland Forest or Shrubland Forest or Shrubland
L1: Protected or
IRA Status
Not Protected Not Protected Not Protected
L2: Slope <40% <60% <60%
L3: Distance to
Nearest Road
<1,000ft <2,000ft <2,000ft
L4: Riparian
Buffer
>100ft >100ft >50ft
L5: Administrative
Designation
No Designation No Designation Any Designation

The figure below illustrates the workflow used to quantify the amount of land available for mechanical forest health and risk reduction fuel treatments by considering layered operational constraints. This example illustrates scenario 2 for a 10,252 hectare (25,000 acre) fireshed project area near Rustic, Colorado in the Colorado Front Range priority landscape.

4.2 Data Preparation

4.2.1 Load WCS landscapes

WCS priority landscape spatial data from the USFS Geospatial Data Discovery site (accessed 2023-09-13).

# load data
wf_landscapes <- sf::read_sf("../data/Wildfire_Crisis_Strategy_Landscapes/Wildfire_Crisis_Strategy_Landscapes_(Feature_Layer).shp") %>% 
  dplyr::rename_with(tolower) %>% 
  sf::st_make_valid() %>% 
  dplyr::mutate(
    hectares = (as.numeric(sf::st_area(geometry))/10000) 
    , Mil.Hectares = hectares %>% 
      scales::comma(suffix = " M", scale = 1e-6, accuracy = .01)
    , acres = (as.numeric(sf::st_area(geometry))/4046.85642) 
    , Mil.Acres = acres %>% 
      scales::comma(suffix = " M", scale = 1e-6, accuracy = .01)
  ) %>% 
  dplyr::left_join(
    data.frame(state = datasets::state.name, state_abb = datasets::state.abb)
    , by = join_by("state")
  ) %>% 
  dplyr::mutate(
    area_name = paste0(
      ifelse(is.na(state_abb) | state_abb == "", "UNK", state_abb)
      , ": "
      , name
    )
  ) %>% 
  # manual label placement
  dplyr::left_join(
    readr::read_csv("../data/wildfirepriority_WCS202308/area_label_placement.csv")
    , by = dplyr::join_by("area_name")
  )
#rename sf geom column
  names(wf_landscapes)[names(wf_landscapes)==tolower(attr(wf_landscapes, "sf_column"))] = "geometry"
  sf::st_geometry(wf_landscapes) = "geometry"
# set crs
  transform_crs <- sf::st_crs(wf_landscapes)
# view
  wf_landscapes %>% dplyr::glimpse()

4.2.2 Load landscape-level constraint data

Landscape-level constraint analysis data (i.e., row unique by WCS landscape) was created via this Google Earth Engine script.

# load landscape constraint scenario data
constrained_by_scnro_ls <-
  list.files("../data/wildfirepriority_WCS202308/landscapes/",pattern = "\\.csv$") %>%
  purrr::keep(stringr::str_starts(.,"wfpriority_all_sc")) %>% 
  purrr::map(function(x){
    readr::read_csv(
      paste0("../data/wildfirepriority_WCS202308/landscapes/",x)
      , name_repair = "universal"
      , col_types = cols(.default = "c")
    ) %>% 
    dplyr::rename_with(tolower) %>% 
    dplyr::rename_with(make.names) %>% 
    dplyr::select(state,name,tidyselect::ends_with("_m2"),tidyselect::starts_with("pct_")) %>% 
    dplyr::mutate(scenario_id = stringr::word(x,3,sep="_") %>% readr::parse_number()) %>% 
    dplyr::relocate(scenario_id)
  }) %>%
  dplyr::bind_rows() %>%
  dplyr::inner_join(
    wf_landscapes %>%
      sf::st_drop_geometry() %>%
      dplyr::select(state,name,area_name)
    , by = dplyr::join_by(state,name)
  ) %>%
  dplyr::relocate(area_name,.after = "scenario_id") %>%
  dplyr::group_by(area_name,scenario_id) %>%
  dplyr::filter(dplyr::row_number()==1) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(
    dplyr::across(
      tidyselect::ends_with("_m2")
      , ~ as.numeric(.x) / 10000
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ as.numeric(.x)
    )
  ) %>%
  dplyr::rename_with(
    ~ gsub("_m2", "_ha", .x)
    , tidyselect::ends_with("_m2")
  ) %>%
  # calculate pct reduction
  dplyr::mutate(
    pct_rdctn1_protected = -1*(1 - pct_rmn1_protected)
    , pct_rdctn2_slope = -1*(pct_rmn1_protected - pct_rmn2_slope)
    , pct_rdctn3_roads = -1*(pct_rmn2_slope - pct_rmn3_roads)
    , pct_rdctn4_riparian = -1*(pct_rmn3_roads - pct_rmn4_riparian)
    , pct_rdctn5_administrative = -1*(pct_rmn4_riparian - pct_rmn5_administrative)
    , pct_rdctn_total = -1*(1 - pct_rmn5_administrative)
    , pct_covertype_area = covertype_area_ha/feature_area_ha
    # scenario description
    , scenario_lab = factor(
      scenario_id
      , levels = 1:3
      , labels = paste0("Scenario ", 1:3)
      , ordered = T
    ) %>% forcats::fct_rev()
    , scenario_desc = factor(
      scenario_id
      , levels = 1:3
      , labels = c("Scenario 1\n(status quo)", "Scenario 2\n(slopes+roads)", "Scenario 3\n(slopes+roads+admin)")
      , ordered = T
    )
  )

4.2.3 Load fireshed project area data

The fireshed registry spatial data was obtained from the USFS Geospatial Data Discovery tool for firesheds and project areas (accessed 2023-05-03).

See this section for exploration of fireshed and project area spatial data

### fireshed spatial data
fireshed <- sf::st_read("../data/firesheds/Fireshed_Registry3A_Fireshed/Fireshed_Registry%3A_Fireshed_(Feature_Layer).shp") %>%
  sf::st_transform(transform_crs) %>% 
  setNames(c(
      "shape_id"
      , "area_ha"
      , "fireshed_id"
      , "fireshed_name"
      , "fireshed_code"
      , "fireshed_state"
      , "nopas"
      , "objectid"
      , "fshed_id"
      , "exp_total"
      , "exp_usfs"
      , "exp_nonfs"
      , "exp_usfs_protected"
      , "exp_nonfs_protected"
      , "exp_usfs_managed"
      , "exp_nonfs_managed"
      , "exp_usfs_forest"
      , "exp_nonfs_forest"
      , "exp_usfs_nonforest"
      , "exp_nonfs_nonforest"
      , "exp_usfs_conifer"
      , "exp_nonfs_conifer"
      , "exp_usfs_managedforest"
      , "exp_nonfs_managedforest"
      , "exp_usfs_managedconifer"
      , "exp_nonfs_managedconifer"
      , "exp_nonfs_nonconifer_hihaz"
      , "dist_vs"
      , "crisis_strategy"
      , "key_preformance_indicator"
      , "national_usfs_rank"
      , "national_all_land_rank"
      , "regional_usfs_rank"
      , "regional_all_land_rank"
      , "start_date"
      , "end_date"
      , "geometry"
  )) %>% 
  dplyr::mutate(
    exposure_pct_rank = dplyr::percent_rank(exp_total)
    , exposure_pct_rank_grp = dplyr::case_when(
      exposure_pct_rank >= 1-0.01 ~ "Top 1%"
      , exposure_pct_rank >= 1-0.05 ~ "Top 5%"
      , exposure_pct_rank >= 1-0.10 ~ "Top 10%"
      , exposure_pct_rank >= 1-0.25 ~ "Top 25%"
      , TRUE ~ "Bottom 75%"
    ) %>% 
    factor(
      levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
      , ordered = T
    )
    # there is also a national_all_land_rank column
    , ntllandrank_pct_rank = dplyr::percent_rank(-national_all_land_rank)
    , ntllandrank_pct_rank_grp = dplyr::case_when(
        ntllandrank_pct_rank >= 1-0.01 ~ "Top 1%"
        , ntllandrank_pct_rank >= 1-0.05 ~ "Top 5%"
        , ntllandrank_pct_rank >= 1-0.10 ~ "Top 10%"
        , ntllandrank_pct_rank >= 1-0.25 ~ "Top 25%"
        , TRUE ~ "Bottom 75%"
      ) %>% 
      factor(
        levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
        , ordered = T
      )
    , crisis_strategy = ifelse(is.na(crisis_strategy),"Not High Risk",crisis_strategy) %>% 
      as.factor() %>% 
      forcats::fct_shift()
  )
  #rename sf geom column
    names(fireshed)[names(fireshed)==tolower(attr(fireshed, "sf_column"))] = "geometry"
    sf::st_geometry(fireshed) = "geometry"
    # calculate area
    fireshed <- fireshed %>% 
      dplyr::mutate(
        fireshed_area_ha = as.numeric(sf::st_area(geometry))/10000
        , fireshed_area_acres = (fireshed_area_ha*10000)/4046.85642
      )
## fireshed_proj_area spatial data
fireshed_proj_area <- sf::st_read("../data/firesheds/Fireshed_Registry3A_Project_Area/Fireshed_Registry%3A_Project_Area_(Feature_Layer).shp") %>%
  sf::st_transform(transform_crs) %>% 
  setNames(c(
      "shape_id"
      , "fireshed_id"
      , "pa_id"
      , "pa_area_ha"
      , "objectid"
      , "pa_id2"
      , "fshed_id"
      , "exp_total"
      , "exp_usfs"
      , "exp_nonfs"
      , "exp_usfs_protected"
      , "exp_nonfs_protected"
      , "exp_usfs_managed"
      , "exp_nonfs_managed"
      , "exp_usfs_forest"
      , "exp_nonfs_forest"
      , "exp_usfs_nonforest"
      , "exp_nonfs_nonforest"
      , "exp_usfs_conifer"
      , "exp_nonfs_conifer"
      , "exp_usfs_managedforest"
      , "exp_nonfs_managedforest"
      , "exp_usfs_managedconifer"
      , "exp_nonfs_managedconifer"
      , "exp_nonfs_nonconifer_hihaz"
      , "dist_vs"
      , "pctrecentlydisturbed"
      , "start_date"
      , "end_date"
      , "geometry"
  )) %>% 
  dplyr::mutate(
    exposure_pct_rank = dplyr::percent_rank(exp_total)
    , exposure_pct_rank_grp = dplyr::case_when(
      exposure_pct_rank >= 1-0.01 ~ "Top 1%"
      , exposure_pct_rank >= 1-0.05 ~ "Top 5%"
      , exposure_pct_rank >= 1-0.10 ~ "Top 10%"
      , exposure_pct_rank >= 1-0.25 ~ "Top 25%"
      , TRUE ~ "Bottom 75%"
    ) %>% 
    factor(
      levels = c("Top 1%","Top 5%","Top 10%","Top 25%","Bottom 75%")
      , ordered = T
    )
  )
  #rename sf geom column
    names(fireshed_proj_area)[names(fireshed_proj_area)==tolower(attr(fireshed_proj_area, "sf_column"))] = "geometry"
    sf::st_geometry(fireshed_proj_area) = "geometry"
    # calculate area
    fireshed_proj_area <- fireshed_proj_area %>% 
      dplyr::mutate(
        pa_area_ha = as.numeric(sf::st_area(geometry))/10000
        , pa_area_acres = (pa_area_ha*10000)/4046.85642
      ) %>% 
      # JOIN WITH FIRESHED DATA
      dplyr::inner_join(
        fireshed %>%
          sf::st_drop_geometry() %>%
          dplyr::select(fireshed_id, crisis_strategy, exp_total
                        , exposure_pct_rank, exposure_pct_rank_grp
          ) %>% 
          dplyr::rename(exposure_total=exp_total) %>% 
          dplyr::rename_with(
            ~ paste0("fireshed_",.x)
            , -c(fireshed_id)
          )
        , by = dplyr::join_by(fireshed_id)
      ) %>%
      dplyr::select(pa_id,pa_area_ha
                    ,exp_total,exposure_pct_rank,exposure_pct_rank_grp
                    , tidyselect::starts_with("fireshed_")
      ) %>% 
      dplyr::rename(exposure_total=exp_total) %>%
      dplyr::rename_with(
        ~ paste0("pa_", .x, recycle0 = TRUE)
        , tidyselect::starts_with("exp")
      )

Fireshed project area constraint analysis data (i.e., row unique by WCS landscape + f.p.a.) was created via this Google Earth Engine script

# load fireshed constraint scenario data
constrained_by_scnro_ls_pa <-
  list.files("../data/wildfirepriority_WCS202308/fireshed/",pattern = "\\.csv$") %>%
  purrr::keep(stringr::str_starts(.,"wfpriority_fireshed_sc")) %>% 
  purrr::map(function(x){
    readr::read_csv(
      paste0("../data/wildfirepriority_WCS202308/fireshed/",x)
      , name_repair = "universal"
      , col_types = cols(.default = "c")
    ) %>% 
    dplyr::rename_with(tolower) %>% 
    dplyr::rename_with(make.names) %>% 
    dplyr::select(state,name,pa_id,tidyselect::ends_with("_m2"),tidyselect::starts_with("pct_")) %>% 
    dplyr::mutate(scenario_id = stringr::word(x,3,sep="_") %>% readr::parse_number()) %>% 
    dplyr::relocate(scenario_id)
  }) %>%
  dplyr::bind_rows() %>%
  dplyr::inner_join(
    wf_landscapes %>%
      sf::st_drop_geometry() %>%
      dplyr::select(state,name,area_name)
    , by = dplyr::join_by(state,name)
  ) %>%
  dplyr::relocate(area_name,.after = "scenario_id") %>%
  dplyr::mutate(
    dplyr::across(
      tidyselect::ends_with("_m2")
      , ~ as.numeric(.x) / 10000
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ as.numeric(.x)
    )
  ) %>%
  dplyr::rename_with(
    ~ gsub("_m2", "_ha", .x)
    , tidyselect::ends_with("_m2")
  ) %>%
  dplyr::group_by(area_name,pa_id,scenario_id) %>%
  dplyr::arrange(desc(pct_pa_intrsct)) %>%
  dplyr::filter(dplyr::row_number()==1) %>%
  dplyr::ungroup() %>%
  dplyr::filter(pct_pa_intrsct>=0.00001) %>%
  # calculate pct reduction
  dplyr::mutate(
    pct_rdctn1_protected = -1*(1 - pct_rmn1_protected)
    , pct_rdctn2_slope = -1*(pct_rmn1_protected - pct_rmn2_slope)
    , pct_rdctn3_roads = -1*(pct_rmn2_slope - pct_rmn3_roads)
    , pct_rdctn4_riparian = -1*(pct_rmn3_roads - pct_rmn4_riparian)
    , pct_rdctn5_administrative = -1*(pct_rmn4_riparian - pct_rmn5_administrative)
    , pct_rdctn_total = -1*(1 - pct_rmn5_administrative)
    , pct_covertype_area = covertype_area_ha/feature_area_ha
    # calculate level of constraint
      , cnstrnt_lvl = dplyr::case_when(
          -pct_rdctn_total > 0.8 ~ 1
          , -pct_rdctn_total >= 0.6 ~ 2
          , -pct_rdctn_total >= 0.0 ~ 3
        )
      , cnstrnt_class = factor(
          cnstrnt_lvl 
          , levels = 1:3
          , labels = c("high constraint", "med. constraint", "low constraint")
          , ordered = T
        ) %>% forcats::fct_rev()
      , rmn_cnstrnt_class = factor(
          cnstrnt_lvl 
          , levels = 1:3
          , labels = c("0–19% treatable", "20–40% treatable", ">40% treatable")
          , ordered = T
        ) %>% forcats::fct_rev()
      # scenario description
      , scenario_lab = factor(
          scenario_id
          , levels = 1:3
          , labels = paste0("Scenario ", 1:3)
          , ordered = T
        ) %>% forcats::fct_rev()
      , scenario_desc = factor(
        scenario_id
        , levels = 1:3
        , labels = c("Scenario 1\n(status quo)", "Scenario 2\n(slopes+roads)", "Scenario 3\n(slopes+roads+admin)")
        , ordered = T
      )
  ) %>% 
  # join with fireshed data
  dplyr::inner_join(
    fireshed_proj_area %>%
      sf::st_drop_geometry() %>%
      dplyr::select(pa_id, tidyselect::starts_with("pa_exposure"), tidyselect::starts_with("fireshed_")) %>%
      dplyr::mutate(pa_id=as.character(pa_id))
    , by = dplyr::join_by(pa_id)
  )

4.3 Geographic Delineation Description

Timeline of the research utilized to motivate and support the creation of the WCS:

Nested spatial framework for firesheds. Each scale has specific functionality in terms of the planning processes. Firesheds are the broad scale unit of prioritization, but project areas within them are also prioritized as part of implementation of treatments. Project areas are roughly the size that national forests use for conducting vegetation and fuel management projects. The relative variation among firesheds compared to variation within them controls the relative emphasis on selecting firesheds versus individual planning areas. Ager et al. 2021

#define basemap
states_temp <- USAboundaries::us_states(
  states = c(
        # usfs region 1-6 states
        "MT","WY","CO","NM","AZ","UT","ID","WA","OR","CA","NV"
        # , "KS","NE","SD","ND"
      )
  ) %>% 
  sf::st_transform(transform_crs)
cities_temp <- USAboundaries::us_cities(
  states = c(
        # usfs region 1-6 states
        "MT","WY","CO","NM","AZ","UT","ID","WA","OR","CA","NV"
        # , "KS","NE","SD","ND"
      )
  ) %>% 
  sf::st_transform(transform_crs) %>% 
  dplyr::filter(
    toupper(city) %in% c(
      "DENVER"
      , "TUSCON"
      , "SALT LAKE CITY"
      , "LAS VEGAS"
      , "SAN DIEGO"
      , "LOS ANGELES"
      , "FRESNO"
      , "SAN FRANCISCO"
      , "SACRAMENTO"
      , "PORTLAND"
      , "SEATTLE"
      , "ALBUQUERQUE"
      , "RENO"
      , "BOISE"
      , "HELENA"
      , "BILLINGS"
    )
    | (toupper(city) == "PHOENIX" & state_abbr == "AZ")
  )
# plot basemap
plt_basemap <-
ggplot() + 
  geom_sf(
    data = states_temp
    , fill = NA
    , color = "gray66"
  ) +
  geom_sf_text(
    data = states_temp
    , mapping = aes(label = stusps)
    , color = "gray66"
    , size = 3
  ) +
  geom_sf(
    data = cities_temp
    , shape = 1
    , size = 1
    , color = "gray44"
  ) +
  geom_sf_text(
    data = cities_temp
    , mapping = aes(label = city)
    , color = "gray44"
    , size = 2
    # , hjust = -0.15
    , hjust = 0.5
    , vjust = -0.25
  ) +
  coord_sf(expand = F) +
  theme_void()

# plot basemap simple
plt_basemap_simple <-
ggplot() + 
  geom_sf(
    data = states_temp
    , fill = NA
    , color = "gray66"
  ) +
  coord_sf(expand = F) +
  theme_void()

4.3.1 Fireshed risk to communities (2021)

ggplot() + 
  geom_sf(
    data = fireshed
    , mapping = aes(fill=ntllandrank_pct_rank_grp)
  ) + 
  geom_sf(
    data = USAboundaries::us_states() %>% 
      dplyr::filter(!stusps %in% c("AK","HI","PR")) %>% 
      sf::st_transform(transform_crs)
    , fill = NA
    , color = "gray88"
  ) +
  scale_fill_manual(values=c("firebrick","orange3","khaki","seagreen3","royalblue")) +
  coord_sf(expand = F) +
  labs(
    fill = "Fireshed exposure ranks"
  ) +
  theme_void() + 
  theme(
    legend.position = "top"
    , legend.direction = "horizontal"
    , legend.title = element_text(size = 8)
    , legend.text = element_text(size = 7)
  )

National map of the 7,688 firesheds created from community wildfire transmission data (Evers et al. 2020). The fireshed boundaries were created with a process that delineates hotspots of fire transmission to buildings in adjacent or nearby communities. See the Methods section for details on delineating firesheds. Figure 2 from Ager et al. 2021 (p. 7)

4.3.2 Wilfire Crisis Strategy 21 Priority Landscapes (2022-2023)

plt_fshed <- 
plt_basemap + 
  geom_sf(
    data = fireshed %>%
      sf::st_filter(
        wf_landscapes %>%
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(fill=crisis_strategy)
    , alpha = 0.8
  ) +
  geom_sf(
    data = wf_landscapes
    , mapping = aes(
      color = paste0(investment, " investment")
      # , linetype = paste0(investment, " investment")
    )
    , fill=NA
    # , color="black" 
    , lwd=0.7
  ) +
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = dplyr::case_when(
          stringr::str_starts(name,"Sierra and Elko") ~ name # stringr::word(name,1,3)
          , T ~ 
            stringr::str_wrap(
              dplyr::case_when(
                stringr::str_count(name, "\\w+")<2
                  ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
                , TRUE ~ stringr::word(name,1,2)
              )
            , 7
          )
        )
        , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , stat = "sf_coordinates"
    , size = 2.5
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("gray11", "navy")) + # midnightblue
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Landscapes"
    # , linetype = "Landscapes"
  ) +
  # ggspatial::annotation_scale(location = "br", style = "ticks", pad_y = unit(0.7, "cm"), width_hint = 0.05) +
  theme(
      legend.position = c(0.87,0.87) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 8)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(lwd = 3, fill = NA)) #, linetype = c(5,1,1)
    , fill = guide_legend(order = 2, override.aes = list(lwd = 0))
  )
# plot
plt_fshed

ggplot2::ggsave(
    filename = paste0("../data/plt_map_ls_fshed.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 8.5
    , units = "in"
    , dpi = "print"
    , bg = "white"
  )

The boundaries of the 21 priority landscapes roughly follow the boundaries of “firesheds” prioritized to reduce wildfire transmission to developed areas (Figure 1). Firesheds are geographic delineations with an average area of approximately 250,000 acres (~100,000 hectares) that were created to organize the landscape into units for managing wildfire risk to communities.

4.3.3 Wilfire Crisis Strategy Treatment Objective

plt_basemap + 
  geom_sf(
    data = fireshed %>%
      dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(fill=crisis_strategy)
    , alpha = 0.8
  ) +
  geom_sf(
    data = wf_landscapes
    , mapping = aes(color = paste0(investment, " investment"))
    , fill=NA
    # , color="black" 
    , lwd=0.7
  )+
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = dplyr::case_when(
          stringr::str_starts(name,"Sierra and Elko") ~ name # stringr::word(name,1,3)
          , T ~ 
            stringr::str_wrap(
              dplyr::case_when(
                stringr::str_count(name, "\\w+")<2
                  ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
                , TRUE ~ stringr::word(name,1,2)
              )
            , 7
          )
        )
        , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , stat = "sf_coordinates"
    , size = 2.5
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("black", "navy")) +
  # ggspatial::annotation_scale(location = "br", style = "ticks", pad_y = unit(0.7, "cm"), width_hint = 0.05) +
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Landscapes"
  ) +
  theme(
      legend.position = c(0.87,0.87) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 8)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(lwd = 3, fill = NA)) #, linetype = c(5,1,1)
    , fill = guide_legend(order = 2, override.aes = list(lwd = 0))
  )

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the Strategy (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

4.3.4 Fireshed Project Areas

name_temp <- "Enchanted Circle"
# pa map
plt_fshed_pa <-
ggplot() + 
  geom_sf(data = 
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(
            # fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
            name == name_temp
          ) %>% 
          dplyr::select(pa_id)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      )
    , aes(fill = fireshed_crisis_strategy, color = "Project Area\nboundary")
    , alpha = 0.8
    , lwd = 1.2
  ) +
  geom_sf(
    data = fireshed %>%
      # dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::filter(name == name_temp) %>% 
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(color = "Fireshed\nboundary")
    , fill = NA
    , lwd = 1.2
    , linetype = "dashed"
  ) +
  geom_sf(
    data = wf_landscapes %>% dplyr::filter(name == name_temp)
    , mapping = aes(color = "WCS Landscape\nboundary")
    , fill=NA
    , lwd = 0.9
  )+
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray88")) +
  scale_color_manual(values = c("gray44","skyblue2", "black")) +
  labs(
    fill = "Fireshed\nHigh-Risk Status"
    , color = "Spatial\nFramework"
    , subtitle = name_temp
  ) +
  coord_sf(expand = F) +
  ggspatial::annotation_scale(location = "br", style = "ticks", pad_x = unit(0.1, "cm"), pad_y = unit(0.25, "cm")) +
  theme_void() +
  theme(
      legend.position = c(1.08,0.79) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 10)
      , plot.subtitle = element_text(face = "bold.italic", size = 10, hjust = 0.5)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(size = 2.7, linewidth = 6, fill = NA, linetype = 1)) #, linetype = c(5,1,1)))
    , fill = guide_legend(order = 2, override.aes = list(size = 8,linewidth = 0, alpha = 1)) # , color = NA
  )
# pa map w/in
plt_fshed_pa_incl <- ggplot() + 
  geom_sf(data = 
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(
            # fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
            name == name_temp
          ) %>%
          dplyr::mutate(
            # is_25in = ifelse(pct_pa_intrsct>=0.25,"≥25%","<25%")
            is_25in = ifelse(pct_pa_intrsct>=0.25,"included","excluded")
          ) %>% 
          dplyr::select(
            pa_id
            , is_25in
          )
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      )
    , aes(fill = is_25in, color = "Project Area\nboundary")
    , alpha = 0.9
    , lwd = 1.2
  ) +
  geom_sf(
    data = fireshed %>%
      # dplyr::filter(crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
      sf::st_filter(
        wf_landscapes %>%
          dplyr::filter(name == name_temp) %>% 
          dplyr::select(area_name)
        , .predicate=st_intersects
      )
    , aes(color = "Fireshed\nboundary")
    , fill = NA
    , lwd = 1.2
    , linetype = "dashed"
  ) +
  geom_sf(
    data = wf_landscapes %>% dplyr::filter(name == name_temp)
    , mapping = aes(color = "WCS Landscape\nboundary")
    , fill=NA
    , lwd = 0.9
  )+
  scale_fill_manual(values = viridis::viridis(n=4, direction = -1, alpha = 0.9)[c(1,3)]) +
  scale_color_manual(values = c("gray44","skyblue2", "black")) +
  labs(
    fill = "Project Area\nInclusions"
    # fill = "Included in this analysis\nif ≥25% area w/in landscape"
    , color = "Spatial\nFramework"
    , subtitle = name_temp
  ) +
  coord_sf(expand = F) +
  ggspatial::annotation_scale(location = "br", style = "ticks", pad_x = unit(0.1, "cm"), pad_y = unit(0.25, "cm")) +
  theme_void() +
  theme(
      legend.position = c(1.08,0.825) # "top"
      # , legend.direction = "horizontal"
      , legend.title = element_text(size = 8, face = "bold")
      , legend.text = element_text(size = 10)
      , plot.subtitle = element_text(face = "bold.italic", size = 10, hjust = 0.5)
    ) +
  guides(
    color = guide_legend(order = 1, override.aes = list(size = 2.7, linewidth = 6, fill = NA, linetype = 1)) #, linetype = c(5,1,1)))
    , fill = guide_legend(order = 2, override.aes = list(size = 8,linewidth = 0, alpha = 1)) # , color = NA
  )

# combine
cowplot::plot_grid(
    plotlist = list(
      NULL
      , plt_fshed_pa + 
        theme(
          plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
          , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        )
      , NULL
      , plt_fshed_pa_incl + 
        theme(
          plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
          , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        )
      , NULL
    )
    , nrow = 1
    , rel_widths = c(0.01,1,0.01,1,0.01)
    , labels = c("","A","","B","")
    , label_size = 17
    , label_fontface = "bold"
    , label_colour = "gray11"
    # , label_y = 0.99
  ) +
  theme(plot.background = element_rect(fill = "white", color = NA))

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_ex_spatial_strctr.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 5.8
    , units = "in"
    , dpi = "print"
  )
  1. Nested within firesheds are project areas of approximately 10,000 ha in size which represent the geographic unit at which vegetation and fuel management projects are planned (Ager et al. 2021).

  2. Fireshed project areas that did not have at least 25% of area within the boundary of the Strategy priority landscape were excluded from this analysis. This cutoff was implemented under the assumption that with less than 25% of area available, treatment alone would not substantially affect wildfire behavior across the fireshed (North et al. 2015). For each fireshed project area included in this analysis, the entire area within the project area boundary was used to calculate the percent of mechanically available acreage including area outside of the priority landscape boundary.

4.4 Landscape-level analysis

# filter fireshed project area
constrained_by_scnro_ls_pa <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(pct_pa_intrsct>=0.25)

4.4.1 Reduction Treatable Area Table

# scnum_temp<-2
  tbl_temp <- constrained_by_scnro_ls %>% 
    # dplyr::filter(scenario_id==scnum_temp) %>% 
    dplyr::mutate(
      dplyr::across(
        tidyselect::ends_with("_ha")
        , ~ scales::comma(.x, suffix = "k", scale = 1e-3, accuracy = 1)
      )
      , dplyr::across(
        tidyselect::starts_with("pct_")
        , ~ scales::percent(.x, accuracy = .1)
      )
    ) %>% 
    dplyr::select(
      area_name
      , scenario_lab
      , covertype_area_ha
      , pct_covertype_area
      , pct_rdctn1_protected
      , pct_rdctn2_slope
      , pct_rdctn3_roads
      , pct_rdctn4_riparian
      , pct_rdctn5_administrative
      , rmn5_administrative_area_ha
      , pct_rmn5_administrative
    ) %>% 
    dplyr::arrange(area_name,desc(scenario_lab))
  # make table 
  kableExtra::kable(
      tbl_temp %>% dplyr::select(-c(area_name))
      # , caption = paste0("<b><font color=navy>Scenario "
      #                    ,scnum_temp
      #                    ,"</font></b><br>percent reduction of different types of constraints on mechanical treatment"
      #         )
      , caption = "<b>Wildfire Crisis Strategy 21 Priority Landscapes</b><br>percent reduction of different types of constraints on mechanical treatment"
      , col.names = c(
        ""
        , "Forest+Shrub\n(ha)", "Forest+Shrub\n(%)"
        , "Protected"
        , "Slope\nSteepness"
        , "Road\nDistance"
        , "Riparian\nBuffer"
        , "Administrative"
        , "Mechanically\nAvailable (ha)"
        , "Mechanically\nAvailable (%)"
      )
      , escape = F
    ) %>%  
    add_header_above(c(" " = 1, "Area of\nPriority Landscape"=2, "Constraint\nLeast Flexible to Most Flexible" = 5, " " = 2)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(tbl_temp$area_name))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::scroll_box(width = "740px")
Table 4.2: Wildfire Crisis Strategy 21 Priority Landscapes
percent reduction of different types of constraints on mechanical treatment
Area of
Priority Landscape
Constraint
Least Flexible to Most Flexible
Forest+Shrub (ha) Forest+Shrub (%) Protected Slope Steepness Road Distance Riparian Buffer Administrative Mechanically Available (ha) Mechanically Available (%)
AZ: 4FRI
Scenario 1 2,128k 90.4% -12.4% -7.5% -26.7% -10.2% -15.5% 589k 27.7%
Scenario 2 2,128k 90.4% -12.4% -2.1% -13.3% -14.0% -19.9% 813k 38.2%
Scenario 3 2,128k 90.4% -12.4% -2.1% -13.3% -9.9% 0.0% 1,325k 62.3%
AZ: Prescott
Scenario 1 275k 87.2% -12.9% -16.1% -25.1% -11.9% -4.2% 82k 29.8%
Scenario 2 275k 87.2% -12.9% -3.3% -13.2% -18.6% -6.8% 124k 45.2%
Scenario 3 275k 87.2% -12.9% -3.3% -13.2% -13.2% 0.0% 158k 57.4%
AZ: San Carlos Apache Tribal Forest Protection
Scenario 1 1,119k 91.7% -9.3% -19.7% -49.8% -5.8% -1.1% 161k 14.4%
Scenario 2 1,119k 91.7% -9.3% -5.0% -46.6% -10.1% -2.2% 300k 26.8%
Scenario 3 1,119k 91.7% -9.3% -5.0% -46.6% -7.3% 0.0% 357k 31.9%
CA: North Yuba
Scenario 1 139k 96.5% -11.6% -33.4% -7.8% -14.7% -2.3% 42k 30.2%
Scenario 2 139k 96.5% -11.6% -13.9% -3.2% -23.0% -3.3% 63k 45.0%
Scenario 3 139k 96.5% -11.6% -13.9% -3.2% -16.7% 0.0% 76k 54.6%
CA: Plumas Community Protection
Scenario 1 101k 93.3% 0.0% -23.2% -19.8% -21.5% -1.9% 34k 33.6%
Scenario 2 101k 93.3% 0.0% -5.2% -8.1% -33.9% -2.9% 50k 49.9%
Scenario 3 101k 93.3% 0.0% -5.2% -8.1% -24.4% 0.0% 63k 62.3%
CA: Southern California Fireshed Risk Reduction Strategy
Scenario 1 1,408k 86.0% -57.0% -21.0% -8.3% -2.3% -1.4% 139k 9.9%
Scenario 2 1,408k 86.0% -57.0% -9.6% -6.5% -4.2% -3.0% 278k 19.7%
Scenario 3 1,408k 86.0% -57.0% -9.6% -6.5% -3.0% 0.0% 336k 23.9%
CA: Stanislaus
Scenario 1 116k 94.3% -2.9% -26.7% -10.6% -24.7% -0.6% 40k 34.5%
Scenario 2 116k 94.3% -2.9% -10.7% -3.9% -35.8% -0.8% 53k 46.0%
Scenario 3 116k 94.3% -2.9% -10.7% -3.9% -26.1% 0.0% 66k 56.5%
CA: Trinity Forest Health and Fire-Resilient Rural Communities
Scenario 1 554k 83.5% -21.4% -37.8% -10.4% -8.7% -9.2% 69k 12.4%
Scenario 2 554k 83.5% -21.4% -12.4% -7.8% -17.3% -18.5% 126k 22.7%
Scenario 3 554k 83.5% -21.4% -12.4% -7.8% -12.5% 0.0% 254k 45.9%
CO: Colorado Front Range
Scenario 1 1,052k 72.6% -20.4% -19.9% -20.0% -11.3% -1.6% 281k 26.7%
Scenario 2 1,052k 72.6% -20.4% -5.6% -12.0% -17.3% -3.2% 436k 41.4%
Scenario 3 1,052k 72.6% -20.4% -5.6% -12.0% -12.3% 0.0% 522k 49.6%
ID: Nez Perce-Clearwater-Lower Salmon
Scenario 1 700k 89.0% -43.0% -19.7% -9.1% -1.4% -0.3% 185k 26.5%
Scenario 2 700k 89.0% -43.0% -6.7% -5.5% -2.3% -0.7% 293k 41.8%
Scenario 3 700k 89.0% -43.0% -6.7% -5.5% -1.6% 0.0% 302k 43.2%
ID: Southwest Idaho
Scenario 1 611k 87.6% -13.3% -28.2% -10.4% -2.0% -0.3% 279k 45.8%
Scenario 2 611k 87.6% -13.3% -7.0% -6.6% -2.8% -0.8% 425k 69.5%
Scenario 3 611k 87.6% -13.3% -7.0% -6.6% -2.0% 0.0% 435k 71.2%
MT: Kootenai Complex
Scenario 1 372k 91.4% -21.4% -15.0% -7.8% -6.9% -10.6% 142k 38.3%
Scenario 2 372k 91.4% -21.4% -3.6% -2.1% -8.6% -13.8% 188k 50.6%
Scenario 3 372k 91.4% -21.4% -3.6% -2.1% -6.1% 0.0% 249k 66.9%
NM: Enchanted Circle
Scenario 1 506k 85.6% -8.5% -21.5% -33.4% -7.5% -2.2% 136k 26.8%
Scenario 2 506k 85.6% -8.5% -5.2% -25.7% -11.7% -3.7% 229k 45.2%
Scenario 3 506k 85.6% -8.5% -5.2% -25.7% -8.2% 0.0% 265k 52.3%
NV: Sierra and Elko Fronts
Scenario 1 1,001k 73.1% -26.6% -9.8% -26.8% -6.2% -0.9% 297k 29.6%
Scenario 2 1,001k 73.1% -26.6% -2.0% -15.7% -8.9% -1.5% 453k 45.3%
Scenario 3 1,001k 73.1% -26.6% -2.0% -15.7% -6.3% 0.0% 495k 49.4%
OR: Central Oregon
Scenario 1 959k 87.1% -7.0% -2.8% -15.3% -5.4% -10.0% 571k 59.5%
Scenario 2 959k 87.1% -7.0% -0.6% -5.0% -6.4% -11.3% 669k 69.7%
Scenario 3 959k 87.1% -7.0% -0.6% -5.0% -4.5% 0.0% 796k 83.0%
OR: Klamath River Basin
Scenario 1 2,675k 77.0% -19.3% -14.9% -17.7% -5.4% -4.4% 1,023k 38.3%
Scenario 2 2,675k 77.0% -19.3% -4.9% -7.9% -8.3% -7.8% 1,386k 51.8%
Scenario 3 2,675k 77.0% -19.3% -4.9% -7.9% -5.9% 0.0% 1,659k 62.0%
OR: Mount Hood Forest Health and Fire-Resilient Communities
Scenario 1 326k 76.3% -20.1% -14.3% -17.1% -9.0% -15.7% 78k 23.8%
Scenario 2 326k 76.3% -20.1% -4.1% -7.5% -13.2% -20.1% 114k 35.0%
Scenario 3 326k 76.3% -20.1% -4.1% -7.5% -9.5% 0.0% 191k 58.8%
UT: Pine Valley
Scenario 1 153k 94.2% -37.2% -6.9% -21.5% -6.7% 0.0% 43k 27.8%
Scenario 2 153k 94.2% -37.2% -1.4% -10.5% -9.2% 0.0% 64k 41.7%
Scenario 3 153k 94.2% -37.2% -1.4% -10.5% -6.5% 0.0% 68k 44.4%
UT: Wasatch
Scenario 1 381k 89.5% -43.0% -12.1% -14.8% -4.9% -0.3% 95k 25.0%
Scenario 2 381k 89.5% -43.0% -3.7% -7.9% -6.7% -0.5% 146k 38.3%
Scenario 3 381k 89.5% -43.0% -3.7% -7.9% -4.8% 0.0% 155k 40.7%
WA: Central Washington Initiative
Scenario 1 893k 69.8% -25.0% -29.5% -11.5% -8.1% -12.7% 117k 13.1%
Scenario 2 893k 69.8% -25.0% -9.5% -7.0% -13.4% -22.2% 204k 22.9%
Scenario 3 893k 69.8% -25.0% -9.5% -7.0% -9.7% 0.0% 436k 48.8%
WA: Colville Northeast Washington Vision
Scenario 1 656k 92.8% -20.1% -18.5% -16.0% -10.1% -1.2% 223k 34.1%
Scenario 2 656k 92.8% -20.1% -3.8% -7.5% -14.4% -2.0% 342k 52.2%
Scenario 3 656k 92.8% -20.1% -3.8% -7.5% -10.4% 0.0% 382k 58.3%

4.4.2 Change in Available by Scenario

constrained_by_scnro_ls %>% 
  dplyr::mutate(scenario_lab = scenario_lab %>% forcats::fct_rev()) %>% 
ggplot(
    mapping = aes(x=scenario_lab,y=pct_rmn5_administrative,group=area_name)
    ) +
    geom_line(mapping=aes(color = area_name), linewidth = 1.5) +
    geom_label(
      mapping=aes(label = scales::percent(pct_rmn5_administrative, scale = 100, accuracy = 1))
      , color = "black"
      , size = 3
      , label.padding = unit(0.15, "lines")
    ) +
    facet_wrap(facets = vars(area_name)
      , ncol = 3
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
      , scales = "free_x"
    ) +
    scale_color_viridis_d(option = "turbo", alpha = 0.8) +
    scale_y_continuous(expand = expansion(mult = c(0.1, 0.1)), labels = scales::percent_format(), breaks = scales::extended_breaks(6)) +
    # scale_x_discrete(labels = scales::label_wrap(20)) +
    labs(
      x = "Most Restrictive \U2192 Least Restrictive"
      , y = "Percent Treatable Area Remaining"
    ) +
    theme_light() +
    theme(
      legend.position = "none"
      , axis.text.y = element_text(size=7)
      , axis.text.x = element_text(size=7.5)
      , strip.text = element_text(color = "black", face = "bold", size = 10)
    )

4.4.3 Reduction by Constraint

# reshape
constrained_by_scnro_ls_long <- constrained_by_scnro_ls %>%
  dplyr::select(state, name, area_name, tidyselect::starts_with("scenario_"), tidyselect::starts_with("pct_rdctn")) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("pct_rdctn")
    , names_to = "constraint"
    , values_to = "pct_rdctn"
    , names_prefix = "pct_rdctn"
    , values_drop_na = F
  ) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , constraint = factor(
        constraint
        , ordered = TRUE
        , levels = c(
          "protected"
          , "slope"
          , "roads"
          , "riparian"
          , "administrative"
          , "total"   
        )
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
          , "Total"
        )
      ) %>% forcats::fct_rev()
  ) %>% 
  dplyr::left_join(
    constrained_by_scnro_ls %>% dplyr::select(state,name,scenario_id,pct_rdctn_total)
    , by = join_by(state,name,scenario_id)
  )
# color pallete removing dark blue
n_temp <- ((constrained_by_scnro_ls_long$constraint %>% unique() %>% length())-1)*2
plasma_custom <- viridis::plasma(n = n_temp)[seq(2,n_temp,2)]
# scales::show_col(plasma_custom)
# plasma_custom
# plot
ggplot(
  data = constrained_by_scnro_ls_long %>% 
      dplyr::filter(constraint!="Total") %>% 
      dplyr::mutate(
        dplyr::across(
          c(scenario_lab, area_name)
          , ~ forcats::fct_rev(.x)
        )
      )
  , mapping = aes(y = area_name, x = pct_rdctn)
) +
  geom_col(
    mapping = aes(fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_rdctn<0.11*-1,pct_rdctn,NA), accuracy = 1)
      , group = constraint
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_long %>% 
      dplyr::filter(constraint=="Total") %>% 
      dplyr::mutate(
        dplyr::across(
          c(scenario_lab, area_name)
          , ~ forcats::fct_rev(.x)
        )
      )
    , mapping = aes(
      y = area_name, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(cols = vars(scenario_lab)) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.04, .3)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x =  element_blank() #
    , strip.text = element_text(color = "black", face = "bold")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

ggplot2::ggsave(
    filename = paste0("../data/plt_cnstrnt_rdctn_ls.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

Adding in bar for remaining area:

constrained_by_scnro_ls_long %>% 
  dplyr::mutate(
    constraint = forcats::fct_recode(constraint,"Available"="Total")
    , pct_rdctn_hack = ifelse(constraint=="Available", 1+pct_rdctn, pct_rdctn)
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      c(constraint)
      , ~ forcats::fct_rev(.x)
    )
  ) %>% 
# plot
ggplot(
  mapping = aes(
    y = scenario_lab
    , x = pct_rdctn_hack
    , fill = constraint
    , group = constraint
  )
) +
  geom_col(
    width = 0.7, alpha = 0.9
  ) +
  geom_vline(xintercept = 0, color = "gray55", linetype = "solid", linewidth = 0.7) +
  geom_text(
    mapping = aes(
      label = scales::percent(
        ifelse(pct_rdctn<0.145*-1 & constraint!="Available",pct_rdctn*-1,NA)
        , accuracy = 1
      )
      , group = constraint
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 1.8
  ) +
  geom_text(
    mapping = aes(
      label = scales::percent(
        ifelse(constraint=="Available",pct_rdctn_hack,NA)
        , accuracy = 1
      )
      , group = constraint
      , fontface = "bold"
    )
    , color = "black", size = 2.7
    , hjust = -0.1
  ) +
  annotate(
    geom = "text"
    , x = -0.78
    , y = 0
    , label = "constrained"
    , color = cols_cnstrnt_avail[1]
    , fontface = "bold"
    , size = 2.5
  ) +
  annotate(
    geom = "text"
    , x = 0.88
    , y = 0
    , label = "available"
    , color = cols_cnstrnt_avail[2]
    , fontface = "bold"
    , size = 2.5
  ) +
  scale_fill_manual(values = c(rev(plasma_custom), cols_cnstrnt_avail[2])) +
  scale_x_continuous(
    limits = c(-1.00,1.02)
    , labels = scales::percent_format()
  ) +
  scale_y_discrete(expand = expansion(mult = c(0.8,0.3))) +
  facet_wrap(
    facets = vars(area_name)
    , ncol = 3
    , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    , scales = "free_y"
  ) +
  labs(
    fill = ""
    , x = "Percent Available and Constrained Relative to Total Forest and Shrub Area" 
    # , x = "Percent of Forest & Shrub Area in Overall Landscape" 
    , subtitle = "" 
    , y = ""
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.ticks.x = element_blank()
    , panel.grid.major = element_blank()
    , panel.grid.minor = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 9)
  ) +
  guides(fill = guide_legend(nrow = 1))

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_cnstrnt_avail.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 11
    , units = "in"
    , dpi = "print"
  )
constrained_by_scnro_ls_long %>% 
  dplyr::mutate(
    constraint = forcats::fct_recode(constraint,"Available"="Total")
    , pct_rdctn_hack = ifelse(constraint=="Available", 1+pct_rdctn, pct_rdctn)
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      c(constraint)
      , ~ forcats::fct_rev(.x)
    )
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      c(scenario_lab, area_name)
      , ~ forcats::fct_rev(.x)
    )
  ) %>% 
# plot
ggplot(
  mapping = aes(
    y = area_name
    , x = pct_rdctn_hack
    , fill = constraint
    , group = constraint
  )
) +
  geom_col(
    width = 0.7, alpha = 0.9
  ) +
  geom_vline(xintercept = 0, color = "black", linetype = "solid", linewidth = 1.1) +
  geom_text(
    mapping = aes(
      label = scales::percent(
        ifelse(pct_rdctn<0.145*-1 & constraint!="Available",pct_rdctn*-1,NA)
        , accuracy = 1
      )
      , group = constraint
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2
  ) +
  geom_text(
    mapping = aes(
      label = scales::percent(
        ifelse(constraint=="Available",pct_rdctn_hack,NA)
        , accuracy = 1
      )
      , group = constraint
      , fontface = "bold"
    )
    , color = "black", size = 3
    , hjust = -0.1
  ) +
  annotate(
    geom = "text"
    , x = -0.72
    , y = (constrained_by_scnro_ls_long %>% dplyr::pull(area_name) %>% unique() %>% length())*1.035
    , label = "constrained"
    , color = cols_cnstrnt_avail[1]
    , fontface = "bold"
    , size = 3
  ) +
  annotate(
    geom = "text"
    , x = 0.82
    , y = (constrained_by_scnro_ls_long %>% dplyr::pull(area_name) %>% unique() %>% length())*1.035
    , label = "available"
    , color = cols_cnstrnt_avail[2]
    , fontface = "bold"
    , size = 3
  ) +
  scale_fill_manual(values = c(rev(plasma_custom), cols_cnstrnt_avail[2])) +
  scale_x_continuous(
    limits = c(-1.00,1.02)
    , labels = scales::percent_format()
  ) +
  scale_y_discrete(expand = expansion(mult = c(0.02,0.05))) +
  facet_grid(cols = vars(scenario_lab)) +
  labs(
    fill = ""
    , x = "Percent Available and Constrained Relative to Total Forest and Shrub Area" 
    # , x = "Percent of Forest & Shrub Area in Overall Landscape" 
    , subtitle = "" 
    , y = ""
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.text.y = element_text(color = "black",size=10, face = "bold")
    , axis.ticks.x = element_blank()
    , panel.grid.major = element_blank()
    , panel.grid.minor = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 11)
  ) +
  guides(fill = guide_legend(nrow = 1))

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_cnstrnt_avail_wide.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 6
    , units = "in"
    , dpi = "print"
  )

4.5 Fireshed-level analysis

Based on model simulations of how much area generally needs to be treated to influence wildfire behavior, we binned the firesheds into three classes of mechanical constraint:

  • high (81–100% [i.e., only 0–19% is available for mechanical treatment]): fuels treatment would principally need to rely on fire
  • medium (60–80% [i.e., 20–40% is available]): could use a combination of fire and mechanical thinning
  • low (60% [i.e., >40% is available]): could effectively influence wildfire behavior with mechanical treatment alone

4.5.1 Distribution of Fireshed PA Constraint

constrained_by_scnro_ls_pa %>% 
  dplyr::count(area_name,scenario_lab,cnstrnt_class) %>% 
  dplyr::group_by(area_name,scenario_lab) %>% 
  dplyr::mutate(
    pct = n/sum(n)
    , high_pct=max(ifelse(cnstrnt_class=="high constraint",pct,0)) 
    , tot = sum(n)
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  ) %>% 
ggplot() +
  geom_col(
    mapping = aes(x = pct, y = reorder(area_name, desc(area_name)), fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = reorder(area_name, desc(area_name)), group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.12,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  # geom_text(
  #     data = constrained_byftr_huc12_wide %>%
  #       dplyr::group_by(area_name) %>%
  #       dplyr::summarise(n=n(), high_n=sum(ifelse(cnstrnt_class=="high constraint",1,0)))
  #     , mapping = aes(x = -0.05,y=reorder(area_name, -(high_n/n)),label=paste0("n=",scales::comma(n)))
  #     , size = 3
  #     , color = "black"
  #     , hjust = 0.7
  #   ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(scenario_lab)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "Percent of Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.text.y = element_text(color = "black",size=10, face = "bold")
    , axis.ticks.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 11)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

ggplot2::ggsave(
    filename = paste0("../data/plt_fpa_cnstrnt_dist.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 11
    , units = "in"
    , dpi = "print"
  )

4.5.2 Distribution of Fireshed Area

# plot
ggplot(data = constrained_by_scnro_ls_pa %>% dplyr::filter(scenario_id==1)
         , mapping = aes(x = feature_area_ha, y = reorder(area_name,desc(area_name) ))
    ) +
    geom_vline(
      xintercept = 
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(scenario_id==1) %>% 
          dplyr::pull(feature_area_ha) %>% 
          median()
      , linetype="dashed", color="gray66"
    ) +
    geom_violin(aes(fill = state)) + 
    geom_boxplot(width = 0.15) +
    # geom_point(size = 0.7, color = cols_obj_r[2], alpha = 0.2) +
    geom_text(
      data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(scenario_id==1) %>% 
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(n=n(), max_area=max(feature_area_ha))
      , mapping = aes(
          x = (
            constrained_by_scnro_ls_pa %>% 
            dplyr::filter(scenario_id==1) %>% 
            dplyr::pull(feature_area_ha) %>% 
            min()
          )*.9
          ,y=reorder(area_name,desc(area_name) )
          ,label=paste0("n=",scales::comma(n))
        )
      , size = 3
      , color = "black"
      , hjust = 0.7
    ) +
    scale_fill_viridis_d(option = "viridis", alpha = 0.8) +
    scale_x_continuous(labels = scales::comma, breaks = scales::extended_breaks(n=7)) +
    labs(
      fill = ""
      , y = ""
      , x = "Fireshed Area (ha)"
    ) +
    theme_light() +
    theme(
      legend.position = "none" # c(0.9, 0.9)
      , axis.title = element_text(size=9)
      , axis.text.x = element_text(size=7)
    )

4.5.3 Map of Firesheds by Level of Constraint

# polish for mapping
map_firesheds <- fireshed_proj_area %>% 
  dplyr::select(pa_id, geometry) %>% 
  dplyr::mutate(pa_id=as.character(pa_id)) %>% 
  dplyr::inner_join(
    constrained_by_scnro_ls_pa
    , by = dplyr::join_by("pa_id")
    , multiple = "all"
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ scales::percent(as.numeric(.x), accuracy = 0.1)
    )
    , dplyr::across(
      tidyselect::ends_with("_ha")
      , ~ scales::comma(as.numeric(.x), accuracy = 1)
    )
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  ) %>% 
  dplyr::mutate(
    Priority.Area.Name = area_name
    , Level.of.Constraint = cnstrnt_class
    , Pct.Treatable.Class = rmn_cnstrnt_class
    , Fireshed.ProjArea.ID = pa_id
    , Fireshed.ProjArea.Area.Hectares = feature_area_ha
    , Fireshed.WCS.Status = ifelse(is.na(fireshed_crisis_strategy),"Not High Risk",fireshed_crisis_strategy)
    , ForestShrub.Area.Hectares = covertype_area_ha
    , Pct.Rdctn.Protected = pct_rdctn1_protected
    , Pct.Rdctn.Slope = pct_rdctn2_slope
    , Pct.Rdctn.Roads = pct_rdctn3_roads
    , Pct.Rdctn.Riparian = pct_rdctn4_riparian
    , Pct.Rdctn.Administrative = pct_rdctn5_administrative
    , Treatable.Area.Hectares = rmn5_administrative_area_ha
    , Pct.Treatable.Area = pct_rmn5_administrative
  )
# filter for mapview
mapview_fs_temp <- map_firesheds %>% 
  dplyr::filter(
    scenario_id == 2
  )
# ## filter for high priority only
# hi_map_firesheds <- map_firesheds %>% 
#   dplyr::filter(!is.na(fireshed_crisis_strategy))
# filter for mapview
mapview_fs_sc1_temp <- map_firesheds %>% 
  dplyr::filter(
    scenario_id == 1
  )
mapview_fs_sc3_temp <- map_firesheds %>% 
  dplyr::filter(
    scenario_id == 3
  )
mapview_fs_sc3_hi_temp <- map_firesheds %>% 
  dplyr::filter(
    scenario_id == 3
    & fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")
  )
# mapview
mapview::mapviewOptions(homebutton = FALSE, basemaps = c("OpenStreetMap","Esri.WorldImagery"))
mapview::mapview(wf_landscapes
        , color = "black"
        , lwd = 1
        , alpha.regions = 0
        , label = FALSE
        , legend = FALSE
        , popup = FALSE
        , layer.name = "WCS Landscapes"
        , hide = TRUE
  ) +
mapview::mapview(
  mapview_fs_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 2 Constraints"
  , hide = FALSE
  , popup = leafpop::popupTable(
      mapview_fs_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) + 
mapview::mapview(
  mapview_fs_sc1_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 1 Constraints"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc1_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) + 
mapview::mapview(
  mapview_fs_sc3_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 3 Constraints"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc3_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) +
mapview::mapview(
  mapview_fs_sc3_hi_temp
  , zcol = "cnstrnt_class"
  , col.regions = RColorBrewer::brewer.pal(n=3,name="RdYlBu") %>% rev()
  , alpha.regions = 0.6
  , lwd = 0.2
  , label = FALSE
  , legend = FALSE
  , layer.name = "Scenario 3 Constraints (High Only)"
  , hide = TRUE
  , popup = leafpop::popupTable(
      mapview_fs_sc3_hi_temp
      , zcol = c(
        "Priority.Area.Name"
        , "Level.of.Constraint"
        , "Pct.Treatable.Class"
        , "Fireshed.ProjArea.ID"
        , "Fireshed.ProjArea.Area.Hectares"
        , "Fireshed.WCS.Status"
        , "ForestShrub.Area.Hectares"
        , "Pct.Rdctn.Protected"
        , "Pct.Rdctn.Slope"
        , "Pct.Rdctn.Roads"
        , "Pct.Rdctn.Riparian"
        , "Pct.Rdctn.Administrative"
        , "Treatable.Area.Hectares"
        , "Pct.Treatable.Area"
      )
      , row.numbers = FALSE
      , feature.id = FALSE
    )
) 
# **Scenario 2 used for map above**

4.5.4 Map of Fireshed Project Areas by Scenario

plt_fn_temp <- function(scnum=1) {
  plt_basemap_simple +
    geom_sf(data = map_firesheds %>% dplyr::filter(scenario_id==scnum)
            , aes(fill = cnstrnt_class), lwd = 0.05) +
    
  geom_sf(
    data = wf_landscapes
    # , mapping = aes(color = paste0(investment, " investment"))
    , fill=NA
    , color="black"
    , lwd=0.2
  )+
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = dplyr::case_when(
          stringr::str_starts(name,"Sierra and Elko") ~ name # stringr::word(name,1,3)
          , T ~ 
            stringr::str_wrap(
              dplyr::case_when(
                stringr::str_count(name, "\\w+")<2
                  ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
                , TRUE ~ stringr::word(name,1,2)
              )
            , 7
          )
        )
        # , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , color = "black"
    , stat = "sf_coordinates"
    , size = 2
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  # scale_color_manual(values = c("black", cols_obj_r[2])) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  labs(
    fill = ""
    , tag = paste0(
      map_firesheds %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::pull(scenario_lab) %>% unique()
      , " Map"
    )
  ) +
  theme(
    legend.position = c(0.8,0.9) # "top"
    , legend.text = element_text(size = 7)
    , plot.title = element_text(face = "bold", size = 9, margin = margin(0,0,30,0))
    , plot.tag = element_text(face = "bold", size = 10)
    , plot.tag.position = c(0.1,0.98) # "top"
    # , plot.tag.position = c(0.8,0.94) # "top"
  ) + 
  guides(color = "none")
}
# plot
# plt_fn_temp()
plt_fpas_sc_map <- 
  map_firesheds %>% dplyr::pull(scenario_id) %>% unique() %>% 
    purrr::map(plt_fn_temp)
# print and export maps
  1:length(plt_fpas_sc_map) %>% 
    purrr::map(function(x){
    # save export images for publication
    print(plt_fpas_sc_map[x])
        ggplot2::ggsave(
            filename = paste0("../data/","plt_map_pa_cnstrnt_",x,".jpeg")
            # filename = paste0("../data/terrain/","plt_map_patches_fpa_",x,".jpeg")
            , plot = ggplot2::last_plot()
            , width = 8
            , height = 8.5
            , units = "in"
            , dpi = "print"
            , bg = "white"
          )
      }
    )

4.5.5 Map Fireshed Project Areas by Scenario and Landscape

Define function to map each landscape separately with project areas highlighted by level of constraint. This function is used below to create figures with largest patch highlighted.

#########################################################################
#########################################################################
# Make each plot individually by landscape as solution to small multiples
# this block defines function
#########################################################################
##################hack to align plots for ggmap
ggmap_bbox_fn <- function(map, my_crs=3857) {
    if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
    # Extract the bounding box (in lat/lon) from the ggmap to a numeric vector, 
    # and set the names to what sf::st_bbox expects:
    map_bbox <- setNames(unlist(attr(map, "bb")), c("ymin", "xmin", "ymax", "xmax"))
    # Convert the bbox to an sf polygon, transform it to 3857, 
    # and convert back to a bbox (convoluted, but it works)
    bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), my_crs))
    # Overwrite the bbox of the ggmap object with the transformed coordinates 
    attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
    attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
    attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
    attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
    map
}
plt_crs <- 3857
#########################################################################
#########################################################################
#########################################################################
area_nm_list <- sort(unique(constrained_by_scnro_ls$area_name))
# buffer for smaller areas
if_smaller_ha_temp <- wf_landscapes %>% 
  dplyr::pull(hectares) %>% 
  quantile(0.45) %>% 
  purrr::pluck(1)
# create color palette
# function to plot basemap
landscape_basemap_fn <- function(row_n=1) {
  # should zoom in?
  zoom_level <- dplyr::case_when(
    area_nm_list[row_n] %in% c(
      "WA: Colville Northeast Washington Vision"
      , "CA: Southern California Fireshed Risk Reduction Strategy"
      ) ~ 7
    , area_nm_list[row_n] %in% c(
      "UT: Wasatch"
      , "CA: Trinity Forest Health and Fire-Resilient Rural Communities"
      , "ID: Nez Perce-Clearwater-Lower Salmon"
      , "NV: Sierra and Elko Fronts"
      , "OR: Mount Hood Forest Health and Fire-Resilient Communities"
    ) ~ 8
    , area_nm_list[row_n] %in% c(
      "CA: Plumas Community Protection"
      , "CA: North Yuba"
      , "MT: Kootenai Complex"
      , "UT: Pine Valley"
    ) ~ 9
    , 
      (wf_landscapes %>% 
        dplyr::filter(
            area_name == area_nm_list[row_n]
        ) %>% 
        dplyr::pull(hectares) %>% 
        purrr::pluck(1)
        <= if_smaller_ha_temp
      ) ~ 10
    , TRUE ~ 8
  )
  # should buffer extend?
  buffer_box <- dplyr::case_when(
    area_nm_list[row_n] %in% c("WA: Colville Northeast Washington Vision") ~ 58600
    , area_nm_list[row_n] %in% c(
      "CA: Southern California Fireshed Risk Reduction Strategy"
      , "ID: Nez Perce-Clearwater-Lower Salmon"
      , "MT: Kootenai Complex"
    ) ~ 40000
    , area_nm_list[row_n] %in% c(
      "NV: Sierra and Elko Fronts"
    ) ~ 21000
    , area_nm_list[row_n] %in% c(
      "NM: Enchanted Circle"
      , "UT: Pine Valley"
    ) ~ 25000
    , area_nm_list[row_n] %in% c(
      "AZ: San Carlos Apache Tribal Forest Protection"
      , "CA: Plumas Community Protection"
      , "CA: North Yuba"
      , "CA: Trinity Forest Health and Fire-Resilient Rural Communities"
      , "OR: Mount Hood Forest Health and Fire-Resilient Communities"
    ) ~ 35500
    , area_nm_list[row_n] %in% c("UT: Wasatch") ~ 5000
    , 
      (wf_landscapes %>% 
        dplyr::filter(
            area_name == area_nm_list[row_n]
        ) %>% 
        dplyr::pull(hectares) %>% 
        purrr::pluck(1)
        <= if_smaller_ha_temp
      ) ~ 15500
    , TRUE ~ 6000
  )
  
  # bounding box
  bb_temp <- 
    map_firesheds %>% 
    dplyr::filter(
      area_name == area_nm_list[row_n]
    ) %>% 
    sf::st_union() %>% 
    sf::st_transform(crs=5070) %>% 
    sf::st_buffer(as.numeric(buffer_box)) %>% 
    sf::st_transform(crs=4326) %>% # same as get_map return
    sf::st_bbox()
  # set bbox for get call
  bbox_temp <- c(
    bottom = bb_temp[[2]]
    , top = bb_temp[[4]]
    , right = bb_temp[[3]]
    , left = bb_temp[[1]]
  )
  # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  # stamen maps is being depreciated Oct 31, 2023
  # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  # get map
  hey_ggmap <- ggmap::get_stadiamap(
    bbox = bbox_temp
    , zoom = zoom_level
    , maptype = "stamen_terrain" #"toner-hybrid" # "terrain"
    , crop = T
  )

  # ggmap(hey_ggmap)
  # apply align function
    hey_ggmap_aligned <- ggmap_bbox_fn(hey_ggmap, plt_crs) # Use the function
  # plot
  plt <- ggmap(hey_ggmap_aligned) + 
    coord_sf(
      expand = FALSE
    ) +
    labs(
      title = area_nm_list[row_n]
    ) +
    theme_light() +
    theme(
      legend.position = "bottom"
      , legend.direction  = "horizontal"
      , legend.title = element_text(size=7)
      , legend.margin = margin(c(-10,0,0,0))
      , plot.title = element_text(size = 11)
      , strip.text = element_text(color = "black", face = "bold")
      , axis.title = element_blank()
      , axis.text = element_blank()
      , axis.ticks = element_blank()
      , panel.grid = element_blank()
      , plot.margin = margin(0, 0, 0, 0, "cm")
    )
  # return
  return(plt)
}
# landscape_basemap_fn(row_n = 13)
# function to plot scenario constraints + basemap
scenario_cnstrnt_map_fn <- function(row_n = 1) {
  # location ggspatial::annotation_scale 
  my_scale_loc <- dplyr::case_when(
    area_nm_list[row_n] %in% c(
      "CA: Southern California Fireshed Risk Reduction Strategy"
      , "WA: Colville Northeast Washington Vision"
      , "CO: Colorado Front Range"
      , "AZ: 4FRI"
      , "ID: Southwest Idaho"
    ) ~ "bl"
    , TRUE ~ "br"
  )
  # plot
  plt <- landscape_basemap_fn(row_n) + 
      geom_sf(
        data = map_firesheds %>%
          dplyr::filter(
            area_name == area_nm_list[row_n]
          ) %>%
          sf::st_transform(crs=plt_crs)
        , aes(fill = cnstrnt_class), lwd = 0.05
        , inherit.aes = F
        , alpha = 0.6
      ) +
      geom_sf(
        data = wf_landscapes %>% 
          dplyr::filter(
              area_name == area_nm_list[row_n]
            ) %>% 
          sf::st_transform(crs=plt_crs)
        , fill = NA, color = "black", lwd = 0.3
        , inherit.aes = F
      ) +
      facet_grid(
        cols = vars(scenario_lab)
        , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
        # , switch = "both"
      ) +
      scale_fill_manual(values = cols_rdylbu_lt) +
      ggspatial::annotation_scale(
        location = my_scale_loc
        , style = "ticks"
        , pad_x = unit(0.1, "cm")
        , pad_y = unit(0.1, "cm")
      ) +
      labs(
        fill = "Level of\nConstraint"
      ) +
        guides(
        fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
      )
  # return
  return(plt)
}
# scenario_cnstrnt_map_fn(row_n = 15)
# same output with map of each landscape by scenario but without largest patch highlighted
# go through list of areas to make plot
plt_list <- 1:length(area_nm_list) %>%
  purrr::map(scenario_cnstrnt_map_fn)
# scenario_maps
plt_list[]
####
# area_nm_list
# scenario_cnstrnt_map_fn(18)

4.5.6 Patch Analysis

Notes from Collins et al. 2010:

The larger the individual fuel treatment the greater the potential for tree survival, because of reduced edge effect mortality (Ritchie et al. 2007). Larger individual treatments have a greater potential to reduce fire behavior and slow fire spread, which ultimately impacts adjacent untreated stands and should enhance suppression opportunities and increase firefighter safety (Collins et al. 2010). Planning fewer, larger individual treatments across the landscape appears to be a better strategy when human community protection is a primary concern (Schmidt et al. 2008). Finney et al. (2007) note that the relative improvement of optimized treatment placement breaks down when larger proportions of the landscape (40–50%) are excluded from treatment because of land-management constraints.

# create spatial data of interconnected patches by constraint class
# this is to represent the spatial aggregation of fireshed project areas of similar levels of  mechanical constraint
  # aggregation = Tendency of patch or land-cover types to be spatially adjacent or in close proximity
cnstrnt_class_patches <-
  fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::select(pa_id, area_name, scenario_id, scenario_lab, cnstrnt_class)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      ) %>% 
  # union vectors/polygons by constraint class into one multipolygon
    # this dissolves inner boundaries by constraint class
  dplyr::group_by(area_name, scenario_id, scenario_lab, cnstrnt_class) %>% 
  dplyr::summarize(
    geometry = sf::st_union(geometry, by_feature = F)
  ) %>% 
  # separate a multipolygon geometry into several polygons objects after performing a st_union()
  sf::st_cast(to = "MULTIPOLYGON") %>%
  sf::st_cast(to = "POLYGON") %>% 
  # calculate patch area 
  dplyr::ungroup() %>% 
  dplyr::group_by(area_name, scenario_id, scenario_lab, cnstrnt_class) %>% 
  dplyr::mutate(
    patch_area_ha = as.numeric(sf::st_area(geometry))/10000
    , patch_area_rank = dplyr::dense_rank(dplyr::desc(patch_area_ha))
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::arrange(scenario_id, area_name, cnstrnt_class, dplyr::desc(patch_area_ha)) %>% 
  # calculate proportion of total area in each patch
  dplyr::inner_join(
    # total area of fireshed project areas included for analysis
      # after union to remove internal boundaries
    fireshed_proj_area %>% 
      dplyr::mutate(pa_id=as.character(pa_id)) %>% 
      dplyr::inner_join(
        constrained_by_scnro_ls_pa %>% 
          dplyr::filter(scenario_id==1) %>% 
          dplyr::select(pa_id, area_name)
        , by = dplyr::join_by(pa_id)
        , multiple = "all"
      ) %>% 
      dplyr::group_by(area_name) %>% 
      dplyr::summarize(
        geometry = sf::st_union(geometry, by_feature = F)
      ) %>%
      dplyr::mutate(
        total_area_ha = as.numeric(sf::st_area(geometry))/10000
      ) %>% 
      sf::st_drop_geometry()
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::mutate(
    pct_patch_landscape = patch_area_ha / total_area_ha
  )
# function to highlight largest patch with fireshed project areas under
largest_patch_map_fn <- function(nrow = 1){
    plt <- scenario_cnstrnt_map_fn(nrow) +
      geom_sf(
          data = cnstrnt_class_patches %>%
            dplyr::filter(
              area_name == area_nm_list[nrow]
              & patch_area_rank == 1
              & cnstrnt_class != "med. constraint"
            ) %>%
            sf::st_transform(crs=plt_crs)
          , aes(color = cnstrnt_class, fill = cnstrnt_class)
          , lwd = 0.9
          , inherit.aes = F
          , alpha = 0
          # , show.legend = F
        ) +
      scale_color_manual(values = cols_rdylbu_dk, drop = F) +
      labs(
        color = "Level of\nConstraint"
        , subtitle = "largest patch of interconnected fireshed project areas by level of mechanical constraint" # largest patch of spatially adjacent...
      ) +
      theme(
        # plot.subtitle = element_text(size = 8)
        plot.subtitle = element_blank()
      ) +
      guides(
        color = "none"
        , fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9, linewidth = 0))
      )
  # return
  return(plt)
}
# scenario_cnstrnt_map_fn(4)
# largest_patch_map_fn(5)
# plot each landscape
plt_ptch_list <- 1:length(area_nm_list) %>%
  purrr::map(largest_patch_map_fn)
# scenario_maps
plt_ptch_list[]

# save export images for publication
1:length(plt_ptch_list) %>% 
  purrr::map(function(x){
      # get height vs width
      x_temp = abs(plt_ptch_list[[x]]$data$lon[1]-plt_ptch_list[[x]]$data$lon[4])
      y_temp = abs(plt_ptch_list[[x]]$data$lat[1]-plt_ptch_list[[x]]$data$lat[4])
      # save as long or wide
      ggplot2::ggsave(
          filename = paste0("../data/patch_maps/","plt_map_patches_fpa_",x,".jpeg")
          # filename = paste0("../data/terrain/","plt_map_patches_fpa_",x,".jpeg")
          , plot = plt_ptch_list[[x]] # ggplot2::last_plot()
          , width = ifelse(y_temp/x_temp > 1.5, 8.5, 11) # 11
          , height = ifelse(y_temp/x_temp > 1.5, 11, 8.5) # 8.5
          , units = "in"
          , dpi = "print"
          , bg = "white"
        )
    }
  )

4.5.6.1 Example of how patches were done

name_temp <- "WA: Central Washington Initiative"
# plot project areas and landscape
plt_patchex_stp1 <- ggplot() + 
  geom_sf(
    data = 
      fireshed_proj_area %>% 
        dplyr::mutate(pa_id=as.character(pa_id)) %>% 
        dplyr::select(pa_id) %>% 
        dplyr::inner_join(
          constrained_by_scnro_ls_pa %>% 
            dplyr::filter(
              area_name == name_temp
              & scenario_id == 2
            )
          , by = dplyr::join_by(pa_id)
        )
    , aes(fill = cnstrnt_class), lwd = 0.05
    , inherit.aes = F
    , alpha = 0.6
  ) +
  geom_sf(
    data = wf_landscapes %>% 
      dplyr::filter(
          area_name == name_temp
        ) 
    , fill = NA, color = "black", lwd = 0.3
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  labs(
    fill = "Level of\nConstraint"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , legend.margin = margin(c(-10,0,0,0))
    , plot.title = element_text(size = 11)
    , strip.text = element_text(color = "black", face = "bold")
    , axis.title = element_blank()
    , axis.text = element_blank()
    , axis.ticks = element_blank()
    , panel.grid = element_blank()
    , plot.margin = margin(0, 0, 0, 0, "cm")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9, linewidth = 0))
  )
# plt_patchex_stp1
# plot patches and landscape
plt_patchex_stp2 <-
  ggplot() + 
  geom_sf(
      data = cnstrnt_class_patches %>%
        dplyr::filter(
          area_name == name_temp
          & scenario_id == 2
          # & patch_area_rank == 1
          # & cnstrnt_class != "med. constraint"
        )
      , aes(fill = cnstrnt_class) #, color = cnstrnt_class)
      # , lwd = 1.4
      , alpha = 0.6
    ) +
  geom_sf(
    data = wf_landscapes %>% 
      dplyr::filter(
          area_name == name_temp
        ) 
    , fill = NA, color = "black", lwd = 0.3
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  scale_color_manual(values = cols_rdylbu_dk, drop = F) +
  labs(
    fill = "Level of\nConstraint"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , legend.margin = margin(c(-10,0,0,0))
    , plot.title = element_text(size = 11)
    , strip.text = element_text(color = "black", face = "bold")
    , axis.title = element_blank()
    , axis.text = element_blank()
    , axis.ticks = element_blank()
    , panel.grid = element_blank()
    , plot.margin = margin(0, 0, 0, 0, "cm")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9, linewidth = 0))
  )
# plt_patchex_stp2
# selected largest patch
plt_patchex_stp3 <-
  ggplot() + 
  geom_sf(
      data = cnstrnt_class_patches %>%
        dplyr::filter(
          area_name == name_temp
          & scenario_id == 2
          & patch_area_rank == 1
          & cnstrnt_class != "med. constraint"
        )
      , aes(fill = cnstrnt_class, color = cnstrnt_class)
      , lwd = 1.4
      , alpha = 0.6
    ) +
  geom_sf(
    data = wf_landscapes %>% 
      dplyr::filter(
          area_name == name_temp
        ) 
    , fill = NA, color = "black", lwd = 0.3
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  scale_color_manual(values = cols_rdylbu_dk, drop = F) +
  labs(
    fill = "Level of\nConstraint"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , legend.margin = margin(c(-10,0,0,0))
    , plot.title = element_text(size = 11)
    , strip.text = element_text(color = "black", face = "bold")
    , axis.title = element_blank()
    , axis.text = element_blank()
    , axis.ticks = element_blank()
    , panel.grid = element_blank()
    , plot.margin = margin(0, 0, 0, 0, "cm")
  ) +
  guides(
    color = "none"
    , fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9, linewidth = 0))
  )
# plt_patchex_stp3
### combine
cowplot::plot_grid(
  plotlist = list(plt_patchex_stp1, plt_patchex_stp2, plt_patchex_stp3)
  , cols = 3
  , rows = 1
  , align = "h"
)

ggplot2::ggsave(
    filename = paste0("../data/plt_lrgst_patch_example_methods.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 7
    , units = "in"
    , dpi = "print"
  )

Comparison of largest patch of lightly constrained fireshed project areas versus largest patch of highly constrained fireshed project areas.

Planning fewer, larger individual treatments across the landscape appears to be a better strategy when human community protection is a primary concern (Schmidt et al. 2008). Finney et al. (2007) note that the relative improvement of optimized treatment placement breaks down when larger proportions of the landscape (40–50%) are excluded from treatment because of land-management constraints.

# aggregate to constraint class level
agg_cnstrnt_class_patches <- cnstrnt_class_patches %>% 
  sf::st_drop_geometry() %>% 
  dplyr::group_by(area_name, scenario_lab, cnstrnt_class) %>% 
  dplyr::summarize(
    max_patch_area_ha = max(patch_area_ha, na.rm = T)
    , mean_patch_area_ha = mean(patch_area_ha, na.rm = T)
    , max_pct_patch_landscape = max(pct_patch_landscape, na.rm = T)
    , sum_pct_patch_landscape = sum(pct_patch_landscape)
  ) %>% 
  dplyr::ungroup() %>%
  dplyr::mutate(
    max_pct_patch_landscape_hack = dplyr::case_when(
      cnstrnt_class=="high constraint" ~ -max_pct_patch_landscape
      , TRUE ~ max_pct_patch_landscape
    )
  ) %>% 
  dplyr::filter(
    cnstrnt_class!="med. constraint"
  )
# combine with full data to get 0's
agg_cnstrnt_class_patches <- 
  # full data
  dplyr::cross_join(
    constrained_by_scnro_ls_pa %>% dplyr::select(area_name) %>% dplyr::distinct()
    , constrained_by_scnro_ls_pa %>% dplyr::select(scenario_lab) %>% dplyr::distinct()
  ) %>% 
  dplyr::cross_join(
    constrained_by_scnro_ls_pa %>% dplyr::select(cnstrnt_class) %>% dplyr::distinct()
  ) %>% 
  dplyr::filter(
    cnstrnt_class!="med. constraint"
  ) %>% 
  # end full data
  dplyr::left_join(
    agg_cnstrnt_class_patches
    , by = dplyr::join_by(area_name, scenario_lab, cnstrnt_class)
  ) %>% 
  # fill with 0's
  dplyr::mutate(
    dplyr::across(
      .cols = tidyselect::where(is.numeric) & -c(area_name, scenario_lab, cnstrnt_class, max_pct_patch_landscape_hack)
      , .fn = function(x){ifelse(is.na(x),0,x)}
    )
    , max_pct_patch_landscape_hack = ifelse(is.na(max_pct_patch_landscape_hack),-0.001,max_pct_patch_landscape_hack)
  )
# plot
ggplot(
  data = agg_cnstrnt_class_patches
  , mapping = aes(
    y = scenario_lab
    , x = max_pct_patch_landscape_hack
    , fill = cnstrnt_class
    , group = cnstrnt_class
  )
) +
  # geom_line(
  #   mapping = aes(x = forcats::fct_rev(scenario_lab), y = max_pct_patch_landscape, color = forcats::fct_rev(cnstrnt_class), group = forcats::fct_rev(cnstrnt_class))
  #   , linewidth = 1.2
  # ) +
  geom_col(
    width = 0.7, alpha = 0.9
  ) +
  geom_vline(xintercept = 0, color = "gray55", linetype = "solid", linewidth = 0.7) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack<0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 2.7
    , hjust = +1
  ) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack>=0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 2.7
    , hjust = 0
  ) +
  annotate(
    geom = "text"
    , x = -0.73
    , y = 0
    , label = "high constraint"
    , color = cols_rdylbu_lt[3]
    , fontface = "bold"
    , size = 2.5
  ) +
  annotate(
    geom = "text"
    , x = 0.92
    , y = 0
    , label = "low constraint"
    , color = cols_rdylbu_lt[1]
    , fontface = "bold"
    , size = 2.5
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  scale_x_continuous(
    limits = c(-1.05,1.22)
    , labels = scales::percent_format()
    # , expand = expansion(mult = c(0, 0.1))
    # , breaks = scales::extended_breaks(5)
  ) +
  scale_y_discrete(expand = expansion(mult = c(0.8,0.3))) +
  facet_wrap(
    facets = vars(area_name)
    , ncol = 3
    , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    , scales = "free_y"
  ) +
  labs(
    fill = "Level of\nConstraint"
    , subtitle = "" # "percentage of total landscape area comprised by the largest patch" # "Largest Patch % of Landscape Area"
    , x = "Percent of Overall Landscape Area Comprised by the Largest Patch" # "Largest Patch % of Landscape Area"
    , y = ""
  ) +
  theme_light() +
  theme(
    legend.position = "none" # "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.ticks.x = element_blank()
    , panel.grid.major = element_blank()
    , panel.grid.minor = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 9)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_lrgst_patch_index.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 11
    , units = "in"
    , dpi = "print"
  )

plot by scenario wide for presentation

# median for plot
median_temp <- agg_cnstrnt_class_patches %>%
        dplyr::mutate(
          dplyr::across(
            c(scenario_lab, area_name)
            , ~ forcats::fct_rev(.x)
          )
        ) %>% 
        dplyr::group_by(scenario_lab,cnstrnt_class) %>% 
        dplyr::summarise(
          median_pct = median(max_pct_patch_landscape_hack,na.rm=T)
        ) %>% 
        dplyr::ungroup()
# plot
ggplot(
  data = agg_cnstrnt_class_patches %>%
    dplyr::mutate(
      dplyr::across(
        c(scenario_lab, area_name)
        , ~ forcats::fct_rev(.x)
      )
    )
  , mapping = aes(
    y = area_name
    , x = max_pct_patch_landscape_hack
    , fill = cnstrnt_class
    , group = cnstrnt_class
  )
) +
  geom_vline(
    # median line
    data = median_temp
    , mapping = aes(xintercept=median_pct, group = cnstrnt_class, color = cnstrnt_class)
    , linetype = "dashed"
    , lwd = 0.5
  ) +
  geom_text(
    # median line
    data = median_temp
    , mapping = aes(
      x=median_pct
      , y=0
      , group = cnstrnt_class
      , color = cnstrnt_class
      , label = scales::percent(median_pct, accuracy = 1)
      , fontface = "bold"
    )
    , size = 2
    , hjust = 1
  ) +
  geom_col(
    width = 0.7, alpha = 0.9
  ) +
  geom_vline(xintercept = 0, color = "gray55", linetype = "solid", linewidth = 0.7) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack<0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 3
    , hjust = +1
  ) +
  geom_text(
    mapping = aes(
      label = ifelse(
        max_pct_patch_landscape_hack>=0
        , scales::percent(max_pct_patch_landscape, accuracy = 1)
        , ""
      )
      , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
      , fontface = "bold"
    )
    , color = "black", size = 3
    , hjust = 0
  ) +
  annotate(
    geom = "text"
    , x = -0.75
    , y = 0
    , label = "high constraint"
    , color = cols_rdylbu_lt[3]
    , fontface = "bold"
    , size = 2
  ) +
  annotate(
    geom = "text"
    , x = 0.95
    , y = 0
    , label = "low constraint"
    , color = cols_rdylbu_lt[1]
    , fontface = "bold"
    , size = 2
  ) +
  scale_fill_manual(values = cols_rdylbu_lt) +
  scale_color_manual(values = cols_rdylbu_lt) +
  scale_x_continuous(
    limits = c(-1.05,1.22)
    , labels = scales::percent_format()
    # , expand = expansion(mult = c(0, 0.1))
    # , breaks = scales::extended_breaks(5)
  ) +
  scale_y_discrete(expand = expansion(mult = c(0.07,0.02))) +
  facet_grid(
    cols = vars(scenario_lab)
  ) +
  labs(
    fill = "Level of\nConstraint"
    , subtitle = "" # "percentage of total landscape area comprised by the largest patch" # "Largest Patch % of Landscape Area"
    , x = "Percent of Overall Landscape Area Comprised by the Largest Patch" # "Largest Patch % of Landscape Area"
    , y = ""
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.text.y = element_text(color = "black",size=10, face = "bold")
    , axis.ticks.x = element_blank()
    , panel.grid.major = element_blank()
    , panel.grid.minor = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 11)
  )

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_lrgst_patch_index_wide.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 6
    , units = "in"
    , dpi = "print"
  )

Combine map with proportion of overall landscape in largest patch plot for publication

patch_map_bars_combo_fn <- function(area_num_temp = 20) {

  # plot large patch proportion of area horizontal to match map
  plt_patch_prop_temp <- 
    ggplot(
      data = agg_cnstrnt_class_patches %>%
        dplyr::filter(area_name==area_nm_list[area_num_temp]) %>% 
        dplyr::mutate(scenario_lab = forcats::fct_rev(scenario_lab))
      , mapping = aes(
        y = 1
        , x = max_pct_patch_landscape_hack
        , fill = cnstrnt_class
        , group = cnstrnt_class
      )
    ) +
      geom_col(
        width = 0.7, alpha = 0.9
      ) +
      geom_vline(xintercept = 0, color = "gray55", linetype = "solid", linewidth = 0.7) +
      geom_text(
        mapping = aes(
          label = ifelse(
            max_pct_patch_landscape_hack<0
            , scales::percent(max_pct_patch_landscape, accuracy = 1)
            , ""
          )
          , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
          , fontface = "bold"
        )
        , color = "black", size = 4
        , hjust = +1
      ) +
      geom_text(
        mapping = aes(
          label = ifelse(
            max_pct_patch_landscape_hack>=0
            , scales::percent(max_pct_patch_landscape, accuracy = 1)
            , ""
          )
          , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
          , fontface = "bold"
        )
        , color = "black", size = 4
        , hjust = 0
      ) +
      # label ha below pct
      geom_text(
        mapping = aes(
          label = ifelse(
            max_pct_patch_landscape_hack<0
            , paste0(
                scales::comma(max_patch_area_ha,suffix = "k", scale = 1e-3, accuracy = 1)
                , " ha"
            )
            , ""
          )
          , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
          , fontface = "bold"
        )
        , color = "black", size = 2.7
        , hjust = +1
        , vjust = 1.9
      ) +
      geom_text(
        mapping = aes(
          label = ifelse(
            max_pct_patch_landscape_hack>=0
            , paste0(
                scales::comma(max_patch_area_ha,suffix = "k", scale = 1e-3, accuracy = 1)
                , " ha"
            )
            , ""
          )
          , x = max_pct_patch_landscape_hack #+ .07*sign(max_pct_patch_landscape_hack)
          , fontface = "bold"
        )
        , color = "black", size = 2.7
        , hjust = 0
        , vjust = 1.9
      ) +
      annotate(
        geom = "text"
        , x = -0.85
        , y = 0.0
        , label = "high constraint"
        , color = cols_rdylbu_lt[3]
        , fontface = "bold"
        , size = 2.5
      ) +
      annotate(
        geom = "text"
        , x = 0.85
        , y = 0.0
        , label = "low constraint"
        , color = cols_rdylbu_lt[1]
        , fontface = "bold"
        , size = 2.5
      ) +
      scale_fill_manual(values = cols_rdylbu_lt) +
      scale_x_continuous(
        limits = c(-1.22,1.22)
        , labels = scales::percent_format()
        , position = "bottom"
      ) +
      scale_y_discrete() +
      facet_grid(
        cols = vars(scenario_lab)
        , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
        # , switch = "both"
      ) +
      labs(
        fill = "Level of\nConstraint"
        , subtitle = "" # "percentage of total landscape area comprised by the largest patch" # "Largest Patch % of Landscape Area"
        , x = "Percent of Overall Landscape Area Comprised by the Largest Patch" # "Largest Patch % of Landscape Area"
        , y = ""
      ) +
      theme_light() +
      theme(
        legend.position = "none" # "top"
        , legend.direction  = "horizontal"
        , legend.title = element_text(size=7)
        , axis.title.x = element_text(size=8, face = "bold")
        , axis.title.y = element_blank()
        , axis.text.x = element_blank()
        , axis.ticks.x = element_blank()
        , panel.grid.major = element_blank()
        , panel.grid.minor = element_blank()
        , strip.text = element_text(color = "black", face = "bold", size = 8)
        # , strip.text = element_blank()
      ) +
      guides(
        fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
      )
  # combine
  plt <- cowplot::plot_grid(
      plotlist = list(
        plt_ptch_list[[area_num_temp]] + 
          theme(
            legend.position = "top"
            , plot.title = element_blank()
            , plot.margin = margin(0.5,0.1,0,0.1, "cm")
          )
        # , plt_fshed_pa + 
        #   theme(
        #     plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
        #     , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        #   )
        , plt_patch_prop_temp +
          theme(
            plot.margin = margin(0,0.9,0.5,0.9, "cm")
          )
        #######################################
        # , cowplot::plot_grid(
        #   plotlist = list(
        #     NULL
        #     , plt_patch_prop_temp +
        #       theme(
        #         plot.margin = margin(0,0.2,0.5,0.2, "cm")
        #       )
        #     , NULL
        #   )
        #   , nrow = 1
        #   , rel_widths = c(0.1,1,0.1)
        # )
        ######################################
        # , plt_fshed_pa_incl + 
        #   theme(
        #     plot.background = element_rect(color = "gray77", fill="white", linewidth = 1)
        #     , plot.margin = margin(0.5,3,0.5,0.5, "cm")
        #   )
      )
      , ncol = 1
      , axis = "lr"
      , align = "v"
      , rel_heights = c(1, 0.25)
      , labels = c("A","B")
      , label_size = 17
      , label_fontface = "bold"
      , label_colour = "gray11"
      # , label_y = 0.99
    ) +
    theme(plot.background = element_rect(fill = "white", color = NA))
  # return plot
  return(plt)
}
# these are different sizes so manually set size
c(20) %>% 
  purrr::map(function(x){
    # export
    ggplot2::ggsave(
        filename = paste0("../data/plt_ex_patch_map_bars_combo_",x,".jpeg")
        , plot = patch_map_bars_combo_fn(x)
        , width = 7.7
        , height = 7
        , units = "in"
        , dpi = "print"
      )
  })
c(9) %>% 
  purrr::map(function(x){
    # export
    ggplot2::ggsave(
        filename = paste0("../data/plt_ex_patch_map_bars_combo_",x,".jpeg")
        , plot = patch_map_bars_combo_fn(x)
        , width = 7
        , height = 7.3
        , units = "in"
        , dpi = "print"
      )
  })
c(6) %>% 
  purrr::map(function(x){
    # export
    ggplot2::ggsave(
        filename = paste0("../data/plt_ex_patch_map_bars_combo_",x,".jpeg")
        , plot = patch_map_bars_combo_fn(x)
        , width = 9
        , height = 6.4
        , units = "in"
        , dpi = "print"
      )
  })
# c(20,9,6) %>% 
#   purrr::map(function(x){
#     # export
#     ggplot2::ggsave(
#         filename = paste0("../data/plt_ex_patch_map_bars_combo_",x,".jpeg")
#         , plot = patch_map_bars_combo_fn(x)
#         , width = 7.7
#         , height = 7
#         , units = "in"
#         , dpi = "print"
#       )
#   })

4.6 High-risk fireshed project areas

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the plan (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

4.6.1 Project Areas by WCS Status

# aggregate data by area, status
wcs_status_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(scenario_id==1) %>%
  dplyr::group_by(area_name,fireshed_crisis_strategy) %>% 
  dplyr::summarise(
    n = n()
    , feature_area_ha = sum(feature_area_ha)
  ) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(
    tot_n = sum(n)
    , tot_area = sum(feature_area_ha)
    , pct_n = n/tot_n
    , pct_area = feature_area_ha/tot_area
  ) %>% 
  dplyr::ungroup()
# plot
ggplot(data = wcs_status_temp) +
  geom_col(
    mapping = aes(y = reorder(area_name,tot_n), x = n, fill = fireshed_crisis_strategy)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = area_name, x = n, group = fireshed_crisis_strategy
      , label = scales::comma(ifelse(n>20,n,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(scenario_id==1) %>%
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(
          n = n()
          , feature_area_ha = sum(feature_area_ha)
        )
    , mapping = aes(
      y = reorder(area_name,n), x = n
      , label = scales::comma(n, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray")) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.1)),labels = scales::comma_format()) +
  labs(
    fill = "High-Risk\nStatus"
    , y = ""
    , x = "# of Fireshed\nProject Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  )

4.6.2 Area of high-risk project areas

Only including fireshed project areas that have at least 25% of their area within the bounds of the priority landscapes and are in a high-risk fireshed.

See this section for discussion of high-risk firesheds and criteria for including fireshed project areas.

# data filter
wcs_status_temp %>% 
  dplyr::filter(fireshed_crisis_strategy!="Not High Risk") %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(tot_area = sum(feature_area_ha)) %>% 
# plot
ggplot() +
  geom_col(
    mapping = aes(y = reorder(area_name,tot_area), x = feature_area_ha, fill = fireshed_crisis_strategy)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = area_name, x = feature_area_ha, group = fireshed_crisis_strategy
      , label = ifelse(
          feature_area_ha<150000
          ,NA
          ,scales::comma(feature_area_ha,suffix = "k", scale = 1e-3, accuracy = 1)
        )
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_pa %>% 
        dplyr::filter(
          scenario_id==1
          & fireshed_crisis_strategy!="Not High Risk"
        ) %>%
        dplyr::group_by(area_name) %>% 
        dplyr::summarise(
          n = n()
          , feature_area_ha = sum(feature_area_ha)
        )
    , mapping = aes(
      y = reorder(area_name,feature_area_ha), x = feature_area_ha
      , label = scales::comma(feature_area_ha,suffix = "k ha", scale = 1e-3, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  scale_fill_manual(values = c("All-Lands"="firebrick","USFS-Only"="orange3","Not High Risk"="gray")) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.15)),labels = scales::comma_format(suffix = "k ha", scale = 1e-3, accuracy = 1)) +
  labs(
    fill = "High-Risk\nStatus"
    , y = ""
    , x = "Area (ha) of Fireshed\nProject Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
  )

4.6.3 Reduction by Constraint

See this section for reduction by constraint at the landscape-level.

# aggregate by area_name, scenario for hi-risk pa's
constrained_by_scnro_ls_hionly_long <-
  constrained_by_scnro_ls_pa %>%
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha
  ) %>%
  dplyr::group_by(
    scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
  ) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , constraint = factor(
        constraint
        , ordered = TRUE
        , levels = c(
          "protected"
          , "slope"
          , "roads"
          , "riparian"
          , "administrative"
        )
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
        )
      ) %>% forcats::fct_rev()
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  )
  
# plot
ggplot(data = constrained_by_scnro_ls_hionly_long) +
  geom_col(
    mapping = aes(y = scenario_lab, x = pct_rdctn, fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = scenario_lab, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.06*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = constrained_by_scnro_ls_hionly_long %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = scenario_lab, x = -1*(1 - pct_rmn)
      , label = scales::percent(-1*(1 - pct_rmn), accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_wrap(facets = vars(area_name)
      , ncol = 3
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.02, .21)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
    , subtitle = "High-Risk Fireshed Project Areas"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 10) 
  )

The figure above shows the reduction by constraint for only high-risk fireshed project areas with at least 25% of area within the boundary of the Strategy priority landscape.

4.6.4 Reduction by Constraint high-risk vs. total

# combine data and plot
pct_rdctn_temp <- 
constrained_by_scnro_ls_hionly_long %>% 
  dplyr::ungroup() %>% 
  dplyr::select(scenario_lab, area_name, pct_rdctn, constraint_lvl) %>% 
  dplyr::mutate(dta_src = "High-Risk Project Areas") %>% 
  dplyr::bind_rows(
    constrained_by_scnro_ls_long %>% 
      dplyr::ungroup() %>% 
      dplyr::select(scenario_lab, area_name, pct_rdctn, constraint_lvl) %>% 
      dplyr::filter(!is.na(constraint_lvl)) %>% 
      dplyr::mutate(dta_src = "Overall Landscape Area")
  ) %>% 
  dplyr::group_by(scenario_lab, area_name, dta_src) %>% 
  dplyr::mutate(
    pct_rdctn_total = sum(pct_rdctn)
    , constraint = factor(
        constraint_lvl
        , ordered = TRUE
        , levels = 1:5
        , labels = c(
          "Protected"
          , "Slope\nSteepness"
          , "Road\nDistance"
          , "Riparian\nBuffer"
          , "Administrative"
        )
      ) %>% forcats::fct_rev()
    , scenario_lab = scenario_lab %>% forcats::fct_rev()
  )
# plot
ggplot(data = pct_rdctn_temp) +
  geom_col(
    mapping = aes(y = dta_src, x = pct_rdctn, fill = constraint)
    , color = NA, width = 0.7
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = dta_src, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.075*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.3
  ) +
  geom_text(
    data = pct_rdctn_temp %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = dta_src, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(
      cols = vars(scenario_lab)
      , rows = vars(area_name)
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
      , switch = "y"
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0.02, .21)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text.x = element_text(color = "black", face = "bold")
    , strip.text.y.left = element_text(angle = 0)
    , strip.placement = "outside"
    
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

4.7 Objective vs. Treatable

Models suggest that fuels reduction treatments can effectively reduce fire size and severity when 20-40% of the landscape has been treated (e.g., Finney, 2007) and the Strategy (USFS, 2022a) aims to treat this amount of high-risk fireshed area within the 21 priority landscapes. The total area of the 21 priority landscapes is approximately 47 million acres (19 million ha), of which 23.7 million acres (9.6 million ha) are in high-risk firesheds; an objective of treating 20-40% would indicate the need to treat between 4.7 to 9.5 million acres (1.9 to 3.8 million ha).

Using the area of fireshed project areas that have at least 25% of their area within the bounds of the priority landscapes and are in a high-risk fireshed as the treatment area objective (20-40%).

# Area of pa's by level of constraint
constrained_by_scnro_ls_cnstrnt_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::group_by(scenario_id,scenario_lab,state,name,area_name,cnstrnt_class) %>% 
  dplyr::summarise(
    n = n()
    , dplyr::across(
      c(feature_area_ha,covertype_area_ha,rmn5_administrative_area_ha)
      , ~ sum(.x,na.rm = T)
    )
  ) %>% 
  dplyr::group_by(scenario_id,scenario_lab,state,name,area_name) %>% 
  dplyr::mutate(
    # cnstrnt_class = stringr::word(cnstrnt_class) %>% stringr::str_remove_all("\\.")
    total_area_ha = sum(feature_area_ha)
    , total_covertype_area_ha = sum(covertype_area_ha)
    # "the need to treat 20 to 40 percent of the overall fireshed"
    , pct_treatable =  rmn5_administrative_area_ha/total_area_ha
    , objective_lo_treat_area_ha = total_area_ha*.2
    , objective_hi_treat_area_ha = total_area_ha*.4
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::arrange(state,name,scenario_id,cnstrnt_class)

4.7.1 Objective vs. Treatable Table

table_temp <-
  constrained_by_scnro_ls_cnstrnt_temp %>% 
  dplyr::mutate(
    cnstrnt_class = stringr::word(cnstrnt_class) %>% stringr::str_remove_all("\\.")
  ) %>% 
  dplyr::rename(
    remaining_pct=pct_treatable
    , remaining_ha=rmn5_administrative_area_ha
  ) %>%
  dplyr::select(-c(feature_area_ha, covertype_area_ha)) %>% 
  tidyr::pivot_wider(
    names_from = cnstrnt_class
    , values_from = c(n, remaining_pct, remaining_ha)
    , names_glue = "{cnstrnt_class}_{.value}"
    , values_fill = 0
  ) %>% 
  dplyr::mutate(
    tot_remaining_ha = low_remaining_ha+med_remaining_ha
    ,tot_remaining_pct = low_remaining_pct+med_remaining_pct
  ) %>% 
  dplyr::relocate(tot_remaining_ha,.after = med_remaining_ha) %>% 
  dplyr::inner_join(
    wf_landscapes %>% 
      sf::st_drop_geometry() %>% 
      dplyr::select(area_name,hectares) %>% 
      dplyr::rename(landscape_area_ha=hectares)
    , by = dplyr::join_by(area_name)
  )
# make table 
kableExtra::kable(
    table_temp %>% 
      dplyr::select(
        c(
          scenario_lab
          # , landscape_area_ha
          , total_area_ha
          , total_covertype_area_ha
          # light and mod area
          , low_remaining_ha, med_remaining_ha, tot_remaining_ha
          , low_remaining_pct, med_remaining_pct, tot_remaining_pct
        )
      ) %>% 
      dplyr::mutate(
        dplyr::across(
          tidyselect::ends_with("_ha")
          , ~ scales::comma(.x, suffix = "k", scale = 1e-3, accuracy = 1)
        )
        , dplyr::across(
          tidyselect::ends_with("_pct")
          , ~ scales::percent(.x, accuracy = .1)
        )
      )
  , caption = "Wildfire Crisis Strategy High-Risk Fireshed Project Areas<br>mechanically available for treatment in lightly and moderately constrained project areas (≥25% in landscape)"
  , col.names = c(
    ""
    # , "Total\nLandscape"
    , "Total Area (ha)"
    , "Forest+Shrub Area (ha)"
    , "Lightly Constrained", "Moderately Constrained"
    , "Total"
    , "Lightly Constrained", "Moderately Constrained"
    , "Total"
  
  )
  , escape = F
) %>%  
    add_header_above(c(" " = 1, "High-Risk Fireshed\nProject Areas" = 2, "Treatable High-Risk\nArea (ha)"=3, "Treatable High-Risk\nPercent (%)" = 3)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(table_temp$area_name))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::column_spec(
      9 # 10
      , color = "white"
      , background = 
          table_temp %>% 
            dplyr::mutate(
              is_gt20 = dplyr::case_when(
                low_remaining_pct >= .2 ~ cols_obj_html[1]
                , tot_remaining_pct >= .2 ~ cols_obj_html[2]
                , TRUE ~ cols_obj_html[3]
              )
          ) %>% 
          dplyr::pull(is_gt20)
          
    ) %>% 
    kableExtra::scroll_box(width = "740px")
Table 4.3: Wildfire Crisis Strategy High-Risk Fireshed Project Areas
mechanically available for treatment in lightly and moderately constrained project areas (≥25% in landscape)
High-Risk Fireshed
Project Areas
Treatable High-Risk
Area (ha)
Treatable High-Risk
Percent (%)
Total Area (ha) Forest+Shrub Area (ha) Lightly Constrained Moderately Constrained Total Lightly Constrained Moderately Constrained Total
AZ: 4FRI
Scenario 1 1,579k 1,436k 304k 126k 430k 19.3% 8.0% 27.2%
Scenario 2 1,579k 1,436k 535k 69k 605k 33.9% 4.4% 38.3%
Scenario 3 1,579k 1,436k 913k 50k 963k 57.8% 3.2% 61.0%
AZ: Prescott
Scenario 1 300k 263k 34k 30k 63k 11.3% 9.9% 21.1%
Scenario 2 300k 263k 90k 25k 115k 29.9% 8.5% 38.4%
Scenario 3 300k 263k 134k 15k 150k 44.8% 5.0% 49.8%
AZ: San Carlos Apache Tribal Forest Protection
Scenario 1 250k 230k 4k 15k 19k 1.6% 6.1% 7.7%
Scenario 2 250k 230k 28k 50k 78k 11.3% 19.9% 31.2%
Scenario 3 250k 230k 69k 31k 100k 27.6% 12.2% 39.8%
CA: North Yuba
Scenario 1 206k 197k 21k 36k 57k 10.4% 17.3% 27.6%
Scenario 2 206k 197k 75k 13k 88k 36.5% 6.5% 43.0%
Scenario 3 206k 197k 99k 12k 111k 48.4% 5.7% 54.1%
CA: Plumas Community Protection
Scenario 1 157k 148k 12k 34k 46k 7.8% 21.7% 29.5%
Scenario 2 157k 148k 65k 8k 72k 41.2% 4.9% 46.1%
Scenario 3 157k 148k 91k 0k 91k 57.9% 0.0% 57.9%
CA: Southern California Fireshed Risk Reduction Strategy
Scenario 1 1,053k 854k 13k 92k 105k 1.3% 8.7% 10.0%
Scenario 2 1,053k 854k 162k 87k 249k 15.3% 8.3% 23.7%
Scenario 3 1,053k 854k 241k 71k 311k 22.8% 6.7% 29.6%
CA: Stanislaus
Scenario 1 153k 143k 20k 30k 50k 13.1% 19.8% 32.8%
Scenario 2 153k 143k 59k 10k 70k 38.8% 6.7% 45.5%
Scenario 3 153k 143k 80k 4k 84k 52.4% 2.5% 54.9%
CA: Trinity Forest Health and Fire-Resilient Rural Communities
Scenario 1 565k 488k 3k 32k 36k 0.6% 5.7% 6.3%
Scenario 2 565k 488k 49k 48k 97k 8.6% 8.5% 17.1%
Scenario 3 565k 488k 171k 48k 219k 30.3% 8.5% 38.8%
CO: Colorado Front Range
Scenario 1 1,138k 823k 85k 118k 203k 7.5% 10.4% 17.8%
Scenario 2 1,138k 823k 286k 65k 351k 25.1% 5.7% 30.8%
Scenario 3 1,138k 823k 374k 62k 436k 32.9% 5.4% 38.3%
ID: Nez Perce-Clearwater-Lower Salmon
Scenario 1 520k 444k 77k 49k 126k 14.8% 9.5% 24.2%
Scenario 2 520k 444k 199k 13k 212k 38.3% 2.5% 40.8%
Scenario 3 520k 444k 205k 14k 218k 39.4% 2.6% 42.0%
ID: Southwest Idaho
Scenario 1 734k 632k 228k 41k 268k 31.0% 5.5% 36.5%
Scenario 2 734k 632k 408k 15k 423k 55.6% 2.1% 57.6%
Scenario 3 734k 632k 427k 8k 435k 58.1% 1.1% 59.2%
MT: Kootenai Complex
Scenario 1 322k 294k 98k 12k 111k 30.6% 3.9% 34.4%
Scenario 2 322k 294k 140k 9k 149k 43.5% 2.9% 46.4%
Scenario 3 322k 294k 190k 9k 199k 59.1% 2.7% 61.8%
NV: Sierra and Elko Fronts
Scenario 1 678k 453k 99k 55k 154k 14.6% 8.2% 22.8%
Scenario 2 678k 453k 229k 17k 247k 33.8% 2.6% 36.4%
Scenario 3 678k 453k 250k 20k 270k 36.9% 2.9% 39.9%
NM: Enchanted Circle
Scenario 1 318k 285k 33k 42k 75k 10.4% 13.1% 23.5%
Scenario 2 318k 285k 118k 16k 134k 37.1% 4.9% 42.0%
Scenario 3 318k 285k 133k 13k 146k 41.9% 4.0% 45.9%
OR: Central Oregon
Scenario 1 278k 232k 138k 10k 149k 49.7% 3.6% 53.4%
Scenario 2 278k 232k 169k 0k 169k 60.5% 0.0% 60.5%
Scenario 3 278k 232k 203k 0k 203k 73.1% 0.0% 73.1%
OR: Klamath River Basin
Scenario 1 684k 557k 13k 45k 59k 2.0% 6.6% 8.6%
Scenario 2 684k 557k 137k 36k 173k 20.0% 5.2% 25.3%
Scenario 3 684k 557k 220k 35k 254k 32.1% 5.1% 37.2%
OR: Mount Hood Forest Health and Fire-Resilient Communities
Scenario 1 339k 244k 15k 22k 37k 4.5% 6.5% 11.0%
Scenario 2 339k 244k 46k 21k 67k 13.6% 6.3% 19.9%
Scenario 3 339k 244k 131k 10k 141k 38.6% 2.8% 41.5%
UT: Pine Valley
Scenario 1 187k 173k 7k 32k 39k 3.9% 16.9% 20.8%
Scenario 2 187k 173k 36k 24k 60k 19.1% 13.1% 32.2%
Scenario 3 187k 173k 45k 19k 64k 24.0% 10.3% 34.4%
UT: Wasatch
Scenario 1 156k 126k 21k 11k 32k 13.6% 7.1% 20.7%
Scenario 2 156k 126k 45k 8k 52k 28.6% 5.0% 33.7%
Scenario 3 156k 126k 47k 8k 56k 30.4% 5.3% 35.7%
WA: Central Washington Initiative
Scenario 1 903k 654k 10k 37k 47k 1.1% 4.1% 5.2%
Scenario 2 903k 654k 73k 46k 119k 8.1% 5.1% 13.2%
Scenario 3 903k 654k 285k 43k 328k 31.6% 4.7% 36.3%
WA: Colville Northeast Washington Vision
Scenario 1 272k 244k 55k 28k 83k 20.4% 10.2% 30.6%
Scenario 2 272k 244k 118k 10k 128k 43.3% 3.6% 46.9%
Scenario 3 272k 244k 149k 6k 155k 54.7% 2.2% 56.9%

4.7.2 Reduction in treatable area flow

table_temp %>% 
  dplyr::filter(scenario_id==1) %>% 
  dplyr::select(
    c(
      area_name
      , landscape_area_ha
      , total_area_ha
      , total_covertype_area_ha
      # , tot_remaining_ha
      # , low_remaining_ha
    )
  ) %>%
  tidyr::pivot_longer(
    !area_name
    , names_to = "name"
    , values_to = "hectares"
  ) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(order = dplyr::row_number(),y_val=2) %>% 
  dplyr::ungroup() %>% 
  dplyr::cross_join(data.frame(scenario_id = 1:3)) %>% 
  ####
  dplyr::bind_rows(
    table_temp %>% 
      dplyr::select(
        c(
          area_name
          , scenario_id
          , tot_remaining_ha
          , low_remaining_ha
        )
      ) %>%
      tidyr::pivot_longer(
        !c(area_name,scenario_id)
        , names_to = "name"
        , values_to = "hectares"
      ) %>% 
      dplyr::group_by(area_name,scenario_id) %>% 
      dplyr::mutate(
        order = dplyr::row_number()+3
        , y_val = scenario_id
      ) %>% 
      dplyr::ungroup() %>% 
      dplyr::relocate(scenario_id,.after = dplyr::last_col())
  ) %>% 
  dplyr::arrange(area_name,order,y_val ) %>% 
  dplyr::mutate(
    order_lab = factor(
      order
      , levels = 1:5
      , labels = c(
        "Total\nLandscape"
        , "High-risk\nFireshed PAs"
        , "High-risk\nForest+Shrub"
        , "Lt.+Md. Constrained"
        , "Lightly Constrained"
      )
      , ordered = T
    )
  ) %>% 
  dplyr::inner_join(
    wf_landscapes %>% 
      sf::st_drop_geometry() %>% 
      dplyr::select(area_name, hectares) %>% 
      dplyr::rename(tot_hectares=hectares)
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::inner_join(
    table_temp %>% 
      dplyr::filter(scenario_id==1) %>% 
      dplyr::select(area_name, objective_lo_treat_area_ha)
    , by = dplyr::join_by(area_name)
  ) %>% 
  dplyr::mutate(
    y_val = -y_val
    , sizing = ifelse(hectares>tot_hectares,1,hectares/tot_hectares)
    , coloring = dplyr::case_when(
        order %in% 4:5
          & hectares >= objective_lo_treat_area_ha
          ~ 1
        , order %in% 4:5
          & hectares < objective_lo_treat_area_ha
          ~ 2
        , T ~ 3
    )
  ) %>% 
  # dplyr::filter(stringr::str_starts(area_name,"CA")) %>% 
  # View()
ggplot() +
  geom_line(
    mapping = aes(x = order_lab, y = y_val, group = scenario_id)
    , color = "gray"
  ) +
  geom_point(
    mapping = aes(x = order_lab, y = y_val, size = sizing*100, color = factor(coloring))
    , shape = 15
  ) +
  geom_text(
    mapping = aes(
      x = order_lab, y = y_val, group = scenario_id
      , label = scales::comma(hectares,suffix = "k", scale = 1e-3, accuracy = 1)
    )
    , size = 2.5
    , color = "black"
    , vjust = 2.5
    , angle = 90
  ) +
  geom_text(
    mapping = aes(
      x = order_lab, y = y_val, group = scenario_id
      , label = ifelse(
        order == 4, paste0("Scenario", y_val*-1), NA
      )
      , fontface = "italic"
    )
    , size = 2.5
    , color = "black"
    , hjust = 0.1
    , vjust = -1.4
  ) +
  facet_wrap(facets = vars(area_name)
    , ncol = 3
    , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    , scales = "free_x"
  ) +
  scale_color_manual(values = c(cols_obj_html[1],cols_obj_html[3],cols_obj_html[2])) +
  scale_y_continuous(expand = expansion(c(0.3,0.3))) +
  scale_x_discrete(labels = scales::label_wrap(10)) +
  labs(
    x = ""
    , y = ""
    , subtitle = "<span><span style='color:#008b8b;'><b><i>≥20% treatable</i></b></span> | <span style='color:#cd3700;'><b><i><20% treatable</i></b></span></span>"
  ) +
  theme_light() +
  theme(
    legend.position = "none"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.y = element_blank()
    , axis.ticks.y = element_blank()
    , axis.text = element_text(size = 7)
    , panel.grid = element_blank()
    , plot.subtitle = ggtext::element_markdown(size = 10)
    , strip.text = element_text(color = "black", face = "bold", size = 10)
  )

4.7.3 Objective vs. Treatable Plot

obj_trtbl_temp <-
  constrained_by_scnro_ls_cnstrnt_temp %>% 
  dplyr::filter(cnstrnt_class!="high constraint") %>% 
  dplyr::left_join(
    # aggregate by area_name, scenario_lab
    constrained_by_scnro_ls_cnstrnt_temp %>% 
        dplyr::filter(cnstrnt_class!="high constraint") %>% 
        dplyr::group_by(area_name, scenario_id) %>% 
        dplyr::mutate(low=ifelse(cnstrnt_class=="low constraint",pct_treatable,0)) %>% 
        dplyr::summarise(
          pct_treatable_lowmed = sum(pct_treatable)
          , pct_treatable_low = sum(low)
        ) %>% 
        dplyr::mutate(
          treatment_objective = dplyr::case_when(
              pct_treatable_low >= .2 ~ 1
              , pct_treatable_lowmed >= .2 ~ 2
              , TRUE ~ 3
            ) %>% 
            factor(
              levels = 1:3
              , labels = c(
                  ">20% in\nlow constraint"
                  ,">20% in\nlow+medium constraint"
                  , "<20% in\nlow+medium constraint"
              )
              , ordered = T
            )
        ) %>% 
        dplyr::group_by(scenario_id) %>% 
        dplyr::arrange(scenario_id, treatment_objective
                       , desc(pct_treatable_low), desc(pct_treatable_lowmed)
        ) %>% 
        dplyr::mutate(sorting_rank = dplyr::row_number()) %>% 
        dplyr::ungroup()
      , by = dplyr::join_by(area_name, scenario_id)
  ) %>% 
  dplyr::ungroup() %>% 
  dplyr::mutate(
    dplyr::across(
      c(cnstrnt_class, scenario_lab)
      , ~ forcats::fct_rev(.x)
    )
  ) 
# plot
obj_plt_fn <- function(scnum){
plt <- ggplot(
    data = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum)
    , mapping = aes(
        x = pct_treatable
        , y = reorder(area_name, desc(sorting_rank))
      )
  ) +
  geom_rect(mapping = aes(xmin = 0.2, xmax=0.4, ymin = -Inf, ymax = Inf)
    , fill = "gray77"
    , alpha = 0.6
  ) +
  annotate("text"
     , x = 0.46
     , y = obj_trtbl_temp %>% 
          dplyr::filter(scenario_id == scnum & sorting_rank==max(obj_trtbl_temp$sorting_rank)) %>% 
          dplyr::pull(area_name) %>% 
          unique()
     , label = expression(bold("objective")~bold("20-40%"))
     , size = 3, color = "gray65", parse = T
    ) +
  geom_col(
    mapping = aes(fill = cnstrnt_class)
    , color = NA, width = 0.7
  ) +
  geom_text(
    mapping = aes(
      group = cnstrnt_class
      , label = scales::percent(
          ifelse(pct_treatable<.06,NA,pct_treatable)
          , accuracy = 1
        )
      # , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.2
  ) +
  # total
  geom_text(
    data = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::group_by(area_name, scenario_id) %>% 
        dplyr::filter(dplyr::row_number()==1) %>% 
        dplyr::ungroup()
    , mapping = aes(
      x = pct_treatable_lowmed
      , y = reorder(area_name, desc(sorting_rank))
      , label = scales::percent(pct_treatable_lowmed, accuracy = 1) 
      , color = treatment_objective
      , fontface = "bold"
    )
    , size = 3.5
    , hjust = -0.1
    , show.legend = F
  ) +
  # facet_grid(
  #   rows = vars(scenario_lab)
  #   , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
  #   , scales = "free_y"
  #   # , switch = "y"
  # ) +
  scale_fill_manual(values = RColorBrewer::brewer.pal(n=3,name="RdYlBu")[2:3]) +
  scale_color_manual(values = c(
    ">20% in\nlow constraint"=cols_obj_r[1]
    ,">20% in\nlow+medium constraint"=cols_obj_r[2]
    , "<20% in\nlow+medium constraint"=cols_obj_r[3])
  ) +
  # scale_x_continuous(expand = expansion(mult = c(0, .15)),labels = scales::percent_format()) +
  scale_x_continuous(limits = c(0, .78),labels = scales::percent_format(),position = "top") +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "Percent of Total High-Risk Area"
    , subtitle = obj_trtbl_temp %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::pull(scenario_lab) %>% unique()
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=9, face = "bold")
    , axis.text.x = element_blank()
    , axis.ticks.x = element_blank()
    , axis.text = element_text(size = 8)
    , axis.text.y = element_text(color = "black",size=9, face = "bold")
    , panel.grid.major.x = element_blank()
    , panel.grid.minor.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 10)
    , strip.background = element_rect(fill = "gray88")
    # , strip.placement = "outside"
    , plot.margin = margin(0, 0, 0, 0, "cm")
    , plot.subtitle = element_text(vjust = -4,face="bold")
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )
return(plt)
}
cowplot::plot_grid(
  obj_plt_fn(1) +
    theme(
      legend.position = c(0.86,0.8)
      , legend.direction = "vertical"
    )
  , obj_plt_fn(2) +
    theme(legend.position = "none")
  , obj_plt_fn(3) +
    # scale_x_continuous(position = "top") +
    # labs(x = "Percent of Total High-Risk Area") +
    theme(legend.position = "none", plot.margin = margin(0, 0, 0.1, 0, "cm"))
  , ncol = 1
  , align = 'vh'
)

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_obj_sc.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 11
    , units = "in"
    , dpi = "print"
  )

map landscapes shaded by potential to meet 20-40%

# map landscapes shaded by potential to meet 20-40%
plt_fn_temp <- function(scnum=1) {
  # filter and join data
  sf_dta_temp <- wf_landscapes %>% 
  dplyr::select(area_name) %>% 
    # join with spatial data
    dplyr::inner_join(
      obj_trtbl_temp %>% 
        dplyr::filter(scenario_id==scnum) %>% 
        dplyr::group_by(area_name, scenario_id) %>% 
        dplyr::filter(dplyr::row_number()==1) %>% 
        dplyr::ungroup()
      , by = dplyr::join_by("area_name")
    ) 
  # plot
  plt <- plt_basemap_simple +
  geom_sf(
    data = sf_dta_temp
    , mapping = aes(fill = treatment_objective)
    , color="black"
    , lwd=0.1
    , alpha = 0.9
  ) +
  geom_text(
    data = wf_landscapes %>% 
        sf::st_drop_geometry() %>% 
        sf::st_as_sf(coords = c("lon","lat")) %>% 
        sf::st_set_crs(4326) %>% 
        sf::st_transform(transform_crs)
    , aes(
        label = dplyr::case_when(
          stringr::str_starts(name,"Sierra and Elko") ~ name # stringr::word(name,1,3)
          , T ~ 
            stringr::str_wrap(
              dplyr::case_when(
                stringr::str_count(name, "\\w+")<2
                  ~ stringr::word(name,1,stringr::str_count(name, "\\w+"))
                , TRUE ~ stringr::word(name,1,2)
              )
            , 7
          )
        )
        # , color = paste0(investment, " investment")
        , geometry = geometry
        , fontface = "bold.italic"
      )
    , color = "black"
    , stat = "sf_coordinates"
    , size = 2
    # , seed = 11
    # , min.segment.length = 0
    , show.legend = F
  ) +
  scale_fill_manual(values = c(
    ">20% in\nlow constraint"=cols_obj_r[1]
    ,">20% in\nlow+medium constraint"=cols_obj_r[2]
    , "<20% in\nlow+medium constraint"=cols_obj_r[3])
  ) +
  labs(
    fill = ""
    , tag = paste0(
      map_firesheds %>% dplyr::filter(scenario_id==scnum) %>% 
        dplyr::pull(scenario_lab) %>% unique()
      , " Map"
    )
  ) +
  theme(
    legend.position = c(0.8,0.9) # "top"
    , legend.text = element_text(size = 7)
    , plot.title = element_text(face = "bold", size = 9, margin = margin(0,0,30,0))
    , plot.tag = element_text(face = "bold", size = 10)
    , plot.tag.position = c(0.1,0.98) # "top"
    , legend.spacing.y = unit(0.25, 'cm')
    , plot.margin = margin(0,0,0,0,"cm")
  ) +
  guides(fill = guide_legend(byrow = TRUE))
  # return
  return(plt)
}
# plot
# save export images for publication
obj_trtbl_temp$scenario_id %>%
  unique() %>%
  purrr::map(function(x){
      print(plt_fn_temp(x))
      ggplot2::ggsave(
          filename = paste0("../data/","plt_map_objctv_ls_",x,".jpeg")
          # filename = paste0("../data/terrain/","plt_map_patches_fpa_",x,".jpeg")
          , plot = ggplot2::last_plot()
          , width = 8
          , height = 8.5
          , units = "in"
          , dpi = "print"
          , bg = "white"
        )
    }
  )

make table of landscapes shifting availability toward objective

# kable
tbl_temp <- obj_trtbl_temp %>% 
  dplyr::group_by(area_name, scenario_id) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::group_by(scenario_id, treatment_objective) %>% 
  # dplyr::arrange(desc(pct_treatable_low)) %>% 
  dplyr::arrange(area_name) %>% 
  dplyr::mutate(rnk = dplyr::row_number()) %>% 
  dplyr::ungroup() %>% 
  dplyr::select(scenario_lab, area_name, treatment_objective,rnk) %>% 
  tidyr::pivot_wider(
    names_from = scenario_lab
    , values_from = area_name
  ) %>% 
  dplyr::arrange(treatment_objective, rnk)
# kable table
  options(knitr.kable.NA = "") 
  kableExtra::kable(
    tbl_temp %>% dplyr::select(tidyselect::starts_with("Scenario"))
    , caption = " "
    # , col.names = c("")
  ) %>% 
  kable_classic("striped", full_width=T) %>%
  pack_rows(index = table(forcats::fct_inorder(tbl_temp$treatment_objective)))
Table 4.4:
Scenario 1 Scenario 2 Scenario 3
>20% in
low constraint
ID: Southwest Idaho AZ: 4FRI AZ: 4FRI
MT: Kootenai Complex AZ: Prescott AZ: Prescott
OR: Central Oregon CA: North Yuba AZ: San Carlos Apache Tribal Forest Protection
WA: Colville Northeast Washington Vision CA: Plumas Community Protection CA: North Yuba
CA: Stanislaus CA: Plumas Community Protection
CO: Colorado Front Range CA: Southern California Fireshed Risk Reduction Strategy
ID: Nez Perce-Clearwater-Lower Salmon CA: Stanislaus
ID: Southwest Idaho CA: Trinity Forest Health and Fire-Resilient Rural Communities
MT: Kootenai Complex CO: Colorado Front Range
NM: Enchanted Circle ID: Nez Perce-Clearwater-Lower Salmon
NV: Sierra and Elko Fronts ID: Southwest Idaho
OR: Central Oregon MT: Kootenai Complex
OR: Klamath River Basin NM: Enchanted Circle
UT: Wasatch NV: Sierra and Elko Fronts
WA: Colville Northeast Washington Vision OR: Central Oregon
OR: Klamath River Basin
OR: Mount Hood Forest Health and Fire-Resilient Communities
UT: Pine Valley
UT: Wasatch
WA: Central Washington Initiative
WA: Colville Northeast Washington Vision
>20% in
low+medium constraint
AZ: 4FRI AZ: San Carlos Apache Tribal Forest Protection
AZ: Prescott CA: Southern California Fireshed Risk Reduction Strategy
CA: North Yuba UT: Pine Valley
CA: Plumas Community Protection
CA: Stanislaus
ID: Nez Perce-Clearwater-Lower Salmon
NM: Enchanted Circle
NV: Sierra and Elko Fronts
UT: Pine Valley
UT: Wasatch
<20% in
low+medium constraint
AZ: San Carlos Apache Tribal Forest Protection CA: Trinity Forest Health and Fire-Resilient Rural Communities
CA: Southern California Fireshed Risk Reduction Strategy OR: Mount Hood Forest Health and Fire-Resilient Communities
CA: Trinity Forest Health and Fire-Resilient Rural Communities WA: Central Washington Initiative
CO: Colorado Front Range
OR: Klamath River Basin
OR: Mount Hood Forest Health and Fire-Resilient Communities
WA: Central Washington Initiative

4.8 Totals

##################
# Data Prep
##################
# union landscape and pa level data
# aggregate by scenario for hi-risk pa's
pct_rdctn_pas_temp <- 
  constrained_by_scnro_ls_pa %>%
  dplyr::group_by(pa_id,scenario_id) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::filter(fireshed_crisis_strategy %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha, feature_area_ha
  ) %>%
  dplyr::group_by(scenario_id, scenario_desc, scenario_lab) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  ) %>% 
  dplyr::ungroup()
# overall
pct_rdctn_ls_temp <- 
  constrained_by_scnro_ls %>%
  dplyr::select(scenario_id, scenario_desc, scenario_lab
                , state, name, area_name
                , (tidyselect::starts_with("rmn") & tidyselect::ends_with("area_ha"))
                , covertype_area_ha, feature_area_ha
  ) %>%
  dplyr::group_by(scenario_id, scenario_desc, scenario_lab) %>% 
  dplyr::summarise(
    dplyr::across(
      tidyselect::ends_with("area_ha")
      , sum
    )
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("rmn")
    , names_to = "constraint"
    , values_to = "rmn_area_ha"
    , names_prefix = "rmn"
    , values_drop_na = F
  ) %>%
  dplyr::mutate(constraint = stringr::str_remove(constraint,"_area_ha")) %>% 
  tidyr::separate_wider_delim(constraint, "_", names = c("constraint_lvl", "constraint")) %>% 
  dplyr::mutate(
    constraint_lvl = as.numeric(constraint_lvl)
    , pct_rmn = rmn_area_ha/covertype_area_ha
    , pct_rdctn = dplyr::case_when(
      constraint_lvl == 1 ~ -1*(1-pct_rmn)
      , TRUE ~ -1*(dplyr::lag(pct_rmn)-pct_rmn)
    )
    , pct_rdctn_total = max(-1*(1 - ifelse(constraint_lvl == 5, pct_rmn,0)))
  ) %>% 
  dplyr::ungroup()
# combine
pct_rdctn_all_temp <- 
  pct_rdctn_pas_temp %>% 
  dplyr::mutate(dta_src = "High-Risk Project Areas") %>% 
    dplyr::bind_rows(
      pct_rdctn_ls_temp %>% 
        dplyr::mutate(dta_src = "Overall Landscape Area")
    ) %>% 
    dplyr::mutate(
      constraint_lab = factor(
          constraint_lvl
          , ordered = TRUE
          , 1:5
          , labels = c(
            "Protected"
            , "Slope\nSteepness"
            , "Road\nDistance"
            , "Riparian\nBuffer"
            , "Administrative"
          )
        ) %>% forcats::fct_rev()
    ) %>% 
    dplyr::relocate(dta_src)

4.8.1 Reduction Treatable Area Table

tot_table_temp <- pct_rdctn_all_temp %>% 
  dplyr::mutate(constraint_nm = paste0(constraint_lvl,constraint)) %>% 
  dplyr::select(-c(constraint_lvl,constraint,constraint_lab)) %>% 
  tidyr::pivot_wider(
    names_from = constraint_nm
    , values_from = c(rmn_area_ha, pct_rmn, pct_rdctn)
  ) %>% 
  dplyr::mutate(
    pct_covertype_area = covertype_area_ha/feature_area_ha
  ) %>% 
  dplyr::mutate(
    dplyr::across(
      tidyselect::contains("area_ha")
      , ~ scales::comma(.x,suffix = " M", scale = 1e-6, accuracy = 0.1)
    )
    , dplyr::across(
      tidyselect::starts_with("pct_")
      , ~ scales::percent(.x, accuracy = 0.1)
    )
  ) %>% 
  dplyr::arrange(
    desc(dta_src)
    , scenario_id
  )
# make table 
  kableExtra::kable(
      tot_table_temp %>% 
        dplyr::select(
          scenario_lab
          , covertype_area_ha
          , pct_covertype_area
          , pct_rdctn_1protected
          , pct_rdctn_2slope
          , pct_rdctn_3roads
          , pct_rdctn_4riparian
          , pct_rdctn_5administrative
          , rmn_area_ha_5administrative
          , pct_rmn_5administrative
        )
      , caption = paste0("<b>Wildfire Crisis Strategy Priority Landscapes and High-Risk Fireshed Project Areas</b><br>percent reduction of different types of constraints on mechanical treatment")
      , col.names = c(
        ""
        , "Forest+Shrub\n(ha)", "Forest+Shrub\n%"
        , "Protected"
        , "Slope\nSteepness"
        , "Road\nDistance"
        , "Riparian\nBuffer"
        , "Administrative"
        , "Remaining (ha)"
        , "Remaining (%)"
      )
      , escape = F
    ) %>%  
    add_header_above(c(" " = 1, "Strategy Priority\nArea"=2, "Constraints on\nMechanical Treatment" = 5, " " = 2)) %>%
    kable_classic(full_width=T) %>% 
    pack_rows(index = table(forcats::fct_inorder(tot_table_temp$dta_src))) %>% 
    kableExtra::kable_styling(font_size = 11,fixed_thead = TRUE) %>%  
    kableExtra::scroll_box(width = "740px")
Table 4.5: Wildfire Crisis Strategy Priority Landscapes and High-Risk Fireshed Project Areas
percent reduction of different types of constraints on mechanical treatment
Strategy Priority
Area
Constraints on
Mechanical Treatment
Forest+Shrub (ha) Forest+Shrub % Protected Slope Steepness Road Distance Riparian Buffer Administrative Remaining (ha) Remaining (%)
Overall Landscape Area
Scenario 1 16.1 M 82.9% -22.0% -16.9% -19.8% -7.0% -5.6% 4.6 M 28.7%
Scenario 2 16.1 M 82.9% -22.0% -5.2% -12.0% -10.5% -8.4% 6.8 M 41.9%
Scenario 3 16.1 M 82.9% -22.0% -5.2% -12.0% -7.5% 0.0% 8.6 M 53.3%
High-Risk Project Areas
Scenario 1 8.8 M 82.5% -17.4% -21.0% -18.2% -8.5% -6.3% 2.5 M 28.5%
Scenario 2 8.8 M 82.5% -17.4% -6.4% -10.5% -13.0% -9.5% 3.8 M 43.1%
Scenario 3 8.8 M 82.5% -17.4% -6.4% -10.5% -9.3% 0.0% 5.0 M 56.4%

4.8.2 Reduction by Constraint

# plot
# plot
ggplot(data = pct_rdctn_all_temp) +
  geom_col(
    mapping = aes(y = scenario_lab, x = pct_rdctn, fill = constraint_lab)
    , color = NA, width = 0.6
    , alpha = 0.7
  ) +
  geom_text(
    mapping = aes(
      y = scenario_lab, x = pct_rdctn
      , label = scales::percent(ifelse(pct_rdctn<.06*-1,pct_rdctn,NA), accuracy = 1)
      , fontface = "bold"
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  geom_text(
    data = pct_rdctn_all_temp %>% dplyr::filter(constraint_lvl==5)
    , mapping = aes(
      y = scenario_lab, x = pct_rdctn_total
      , label = scales::percent(pct_rdctn_total, accuracy = 1)
    )
    , color = "black", size = 3.5
    , hjust = -0.1
  ) +
  facet_grid(
      cols = vars(dta_src %>% forcats::fct_rev())
      , labeller = label_wrap_gen(width = 35, multi_line = TRUE)
    ) +
  # scale_fill_viridis_d(option = "plasma") +
  scale_fill_manual(values = plasma_custom) +
  scale_x_reverse(expand = expansion(mult = c(0, .12)),labels = scales::percent_format()) +
  labs(
    fill = ""
    , y = ""
    , x = "Constraint Reduction in Treatable Area (%)"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_blank()
    , strip.text.x = element_text(color = "black", face = "bold")
    , strip.placement = "outside"
    
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

4.8.3 Distribution of Fireshed PA Constraint

cnstrnt_dta_temp <- constrained_by_scnro_ls_pa %>% 
  dplyr::group_by(pa_id,scenario_id) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::filter(fireshed_crisis_strategy  %in% c("USFS-Only", "All-Lands")) %>% 
  dplyr::count(scenario_id, scenario_desc, scenario_lab, cnstrnt_class) %>% 
  dplyr::group_by(scenario_desc, scenario_lab) %>% 
  dplyr::mutate(
    pct = n/sum(n)
    , dta_src = "High-Risk Firesheds"
  ) %>% 
  dplyr::bind_rows(
    constrained_by_scnro_ls_pa %>% 
      dplyr::group_by(pa_id,scenario_id) %>% 
      dplyr::filter(dplyr::row_number()==1) %>% 
      dplyr::ungroup() %>% 
      dplyr::count(scenario_id, scenario_desc, scenario_lab, cnstrnt_class) %>% 
      dplyr::group_by(scenario_desc, scenario_lab) %>% 
      dplyr::mutate(
        pct = n/sum(n)
        , dta_src = "Overall Landscape Area"
      )
  ) %>% 
  dplyr::mutate(
    dta_src = dta_src %>% forcats::fct_rev()
  )
# n for notation
nlab_temp <- 
  cnstrnt_dta_temp %>% 
  dplyr::filter(scenario_id==1) %>% 
  dplyr::group_by(dta_src) %>% 
  dplyr::summarise(n = sum(n)) %>% 
  dplyr::mutate(l = paste0(dta_src, " n = ", scales::comma(n))) %>% 
  dplyr::pull(l) %>% 
  paste(collapse = "\n")
# plot base
plt_fshed_cnstrnt_lvl_temp <- ggplot(data=cnstrnt_dta_temp) +
  geom_col(
    mapping = aes(x = pct, y = scenario_lab, fill=cnstrnt_class)
    ,width = 0.7, alpha=0.8
  ) +
  scale_fill_brewer(type = "div", palette = "RdYlBu", direction = -1) +
  facet_grid(cols = vars(dta_src)) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "Level of\nConstraint"
    , y = ""
    , x = "% of Fireshed\nProject Areas"
    # , caption = nlab_temp
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.title = element_text(size=7)
    , axis.title = element_text(size=9)
    , axis.text.x = element_text(size=7)
    , strip.text = element_text(color = "black", face = "bold")
    , plot.caption = element_text(size = 8)
  ) + 
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )
# plot with large labels
plt_fshed_cnstrnt_lvl_temp + 
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.065,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 4 # 2.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = ifelse(pct>=0.065,
          paste0("n="
            ,scales::comma(n, accuracy = 1)
          )
          , NA)
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.5 # 2.8
    , vjust = 2.3
  ) 

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_prjareas_dist_sc.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

4.8.3.1 Combine map and fireshed plot

# plot small for inset
inset_temp <- plt_fshed_cnstrnt_lvl_temp +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = scales::percent(ifelse(pct>=0.105,pct,NA), accuracy = 1)
        , fontface = "bold"
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3.2 # 2.8
  ) +
  geom_text(
    mapping = aes(x = pct, y = scenario_lab, group=cnstrnt_class
        ,label = ifelse(pct>=0.105,
          paste0("n="
            ,scales::comma(n, accuracy = 1)
          )
          , NA)
      )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 2.1 # 2.8
    , vjust = 2.3
  ) +
    cowplot::theme_minimal_grid(9) +
    theme(
      legend.position = "none"
      , plot.background = element_rect(color = "black", fill="white", linewidth = 2)
      , axis.text.x = element_blank()
      , axis.title.x = element_text(size=7)
      , strip.text = element_text(color = "black", face = "bold", size = 9)
    )

# combine
cowplot::ggdraw(
  plt_fpas_sc_map[[2]] + 
    theme(
      legend.position = c(.4,.98)
      , legend.direction = "horizontal"
      , legend.text = element_text(size = 9)
      , plot.background = element_rect(fill="white", color = NA)
      , plot.margin = margin(0, 5.1, -0.7, 0.1, "cm")
    )
  ) +
  cowplot::draw_plot(inset_temp, x=.5, y=.56, width=.45, height=.4) +
  cowplot::draw_plot_label(
    label = c("A", "B"),
    x = c(0.02, 0.5),
    y = c(1, 0.95),
    size = 17,
    fontface = "bold",
    color = "gray11"
  )

# export
ggplot2::ggsave(
    filename = paste0("../data/plt_prjareas_sc_map_inset.jpeg")
    , plot = ggplot2::last_plot()
    , width = 11
    , height = 8.5
    , units = "in"
    , dpi = "print"
  )

4.9 Landcover Class Area

4.9.1 Overall Landscape

What is the landcover classification of the area within the landscape boundary? This GEE script was used to calculate the area by NLCD landcover class.

# nlcd class groups
nlcd_reclass_df_temp <- data.frame(
    class_id = c(11,12,21,22,23,24,31,41,42,43,51,52,71,72,73,74,81,82,90,95) %>% as.character()
    , class_label = c(
      "Water"
      , "Water"
      , "Developed"
      , "Developed"
      , "Developed"
      , "Developed"
      , "Barren"
      , "Forest"
      , "Forest"
      , "Forest"
      , "Shrubland"
      , "Shrubland"
      , "Herbaceous"
      , "Herbaceous"
      , "Herbaceous"
      , "Herbaceous"
      , "Planted/Cultivated" # double-check name
      , "Planted/Cultivated"
      , "Wetlands"
      , "Wetlands"
      )
  ) %>% 
  dplyr::left_join(
    FedData::nlcd_colors() %>% dplyr::rename(class_id=ID) %>% dplyr::rename_with(tolower) %>% dplyr::mutate_all(as.character)
    , by = dplyr::join_by(class_id)
  )
# load data by wf priority 
wf_nlcd_area <- readr::read_csv("../data/wildfirepriority_WCS202308/wfpriority_nlcd_area.csv") %>% 
  dplyr::rename_with(tolower) %>% 
  dplyr::left_join(
    wf_landscapes %>% 
      dplyr::select(name,state,area_name)
    , by = dplyr::join_by(name,state)
  ) %>% 
  tidyr::pivot_longer(
    cols = tidyselect::starts_with("area_m2_nlcd_cl")
    , names_to = "class_id"
    , values_to = "area_m2"
    , names_prefix = "area_m2_nlcd_cl_"
    , values_drop_na = F
  ) %>% 
  dplyr::left_join(
    nlcd_reclass_df_temp
    , by = dplyr::join_by(class_id)
  ) %>% 
  dplyr::mutate(
    class_label = forcats::fct_relevel(class_label, c("Forest","Shrubland")) %>% forcats::fct_rev()
  )
# get color palette
nlcd_colors <- wf_nlcd_area %>% 
  dplyr::group_by(class_id,color,class_label) %>% 
  dplyr::summarise(sum_area = sum(area_m2,na.rm = T)) %>% 
  dplyr::arrange(class_label,desc(sum_area)) %>% 
  dplyr::group_by(class_label) %>% 
  dplyr::filter(dplyr::row_number()==1) %>% 
  dplyr::ungroup() %>% 
  dplyr::pull(color)
# aggregate and plot
wf_nlcd_area %>% 
  dplyr::group_by(area_name,class_label) %>% 
  dplyr::summarise(sum_area_ha = sum(area_m2,na.rm = F)/10000) %>% 
  dplyr::mutate(sum_area_ha = ifelse(is.na(sum_area_ha),0,sum_area_ha)) %>% 
  dplyr::group_by(area_name) %>% 
  dplyr::mutate(
    pct_area = sum_area_ha/sum(sum_area_ha)
    # , area_name = area_name %>% forcats::fct_rev()
  ) %>%
  dplyr::ungroup() %>% 
ggplot(
    mapping = aes(x = pct_area, y = reorder(area_name, desc(area_name)), fill=class_label)
  ) +
  geom_col(width = 0.7, alpha=0.8) +
  geom_text(
    mapping = aes(
      label = scales::percent(ifelse(pct_area>.03,pct_area,NA), accuracy = 1)
    )
    , position = position_stack(vjust = 0.5)
    , color = "black", size = 3
  ) +
  scale_fill_manual(values = nlcd_colors) +
  scale_x_continuous(labels = scales::percent_format()) +
  labs(
    fill = "NLCD\nCover\nType"
    , y = ""
    , x = "Percent of Overall Landscape Area"
  ) +
  theme_light() +
  theme(
    legend.position = "top"
    , legend.direction  = "horizontal"
    , legend.title = element_text(size=7)
    , axis.title.x = element_text(size=10, face = "bold")
    , axis.title.y = element_blank()
    , axis.text.x = element_blank()
    , axis.text.y = element_text(color = "black",size=10, face = "bold")
    , axis.ticks.x = element_blank()
    , strip.text = element_text(color = "black", face = "bold", size = 11)
  ) +
  guides(
    fill = guide_legend(reverse = T, override.aes = list(alpha = 0.9))
  )

ggplot2::ggsave(
    filename = paste0("../data/plt_nlcd_area_ls.jpeg")
    , plot = ggplot2::last_plot()
    , width = 8.5
    , height = 11
    , units = "in"
    , dpi = "print"
  )